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CHAPTER  1 
INTRODUCTION 


Turbomachinery  components  exhibit  more  diverse  and  complex 
structural  behavior  than  most  classes  of  engineering  structures. 
Stress  analysis  of  rotating  propulsion  system  components  has  been 
a  driving  force  in  the  development  of  many  of  the  most  sophisti¬ 
cated  numerical  methods  in  common  use:  substructuring  and  cyclic 
symmetry  techniques,  large-scale  eigensolution  algorithms,  creep 
and  thermoplasticity  models,  and  modal  and  reduced  basis  methods. 
Even  with  the  powerful  analytical  tools  and  software  which  exist 
today,  the  stress  and  vibration  analysis  of  turbomachine  com¬ 
ponents  is  usually  a  challenging  task.  The  most  important  con¬ 
tributors  to  this  analytical  complexity  are: 

c  . 

yo  intricacy  of  the  structure  geometry  and  properties: 

> -0  nonlinearity  and  its  influence  on  other  responses;  and 

->'□  uncertainty  in  properties,  loading,  and  other  variables.  y 

\  _ _ 

Applied  research  in  finite  element  methods  and  numerical  solution 
algorithms  at  the  present  time  is  concerned,  in  large  measure, 
with  addressing  these  problem  areas. 


This  report  addresses  the  issue  of  uncertainty  in  defining  a 
structural  analysis  model  and  interpreting  the  results.  A  bladed 

disk^4Pigtnre~T)'  is  a  useful  example  of  the  sources  of  uncertainty 

/ 

which  may  exist  for  a  single  model . ~xBlade-to-blade  variations 
may  occur  for  overall  dimensions  or  tnickness  profiles  because  of 
the  manufacturing  processes  involved.!  Material  properties,  even 
within  a  batch  of  material,  change  from  point  to  point.  At  the 
blade  roots,  the  connection  between  blade  and  disk  is  slightly 
different  for  each  blade.  Furthermore,  each  of  these  effects  is 
likely  to  change  as  a  result  of  usage  and  wear.  Finally,  the 
external  forces  acting  on  the  system  include  body  (centripetal) 
forces,  surface  pressures,  and  perhaps  contact  forces  (when  a 
shroud  or  platform  damper  is  present) .  With  the  exception  of  the 
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Figure  1.  Bladed  Disk. 
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centripetal  forces,  these  loads  are  usually  difficult  to  define. 
Pressures,  for  instance,  vary  because  of  flow  paths  established 
by  earlier  stages,  and  the  presence  of  struts  and  other  obstruc¬ 
tions. 

All  of  the  effects  mentioned  above  are  statistical  in 
nature.  It  is  common  practice  to  estimate  them,  conservatively 
if  possible,  for  the  purpose  of  analysis,  and  then  to  perform 
extensive  testing  of  the  finished  system.  However,  certain  of 

these  statistical  effects  are  fundamental  to  the  structural 

.  .  1-5 

behavior  of  interest.  For  example,  the  mistunino  which  re¬ 
sults  from  blade-to-blade  property  variations  in  a  bladed  disk 
system  may  help  to  stabilize  the  flutter  behavior  of  the  system 

due  to  mode  localization  effects6,  but  may  have  a  deleterious 

.  .  7 

effect  on  forced  vibration  amplitudes  . 

This  report  presents  the  development  of  probabilistic  meth¬ 
ods  for  the  analysis  of  turbomachinery  components.  A  consider¬ 
able  body  of  work  exists  in  probabilistic  structural  dynamics, 
and  the  present  study  builds  upon  these  concepts.  We  adopt  a 
middle  ground  in  the  complexity  of  our  statistical  technique,  in 
return  for  ease  in  specifying  the  analytical  problem  and  the 
ability  to  solve  large  problems  at  reasonable  cost.  While  the 
statistical  approach  is  relatively  simple,  it  is  consistent  with 
the  level  of  information  which  is  typically  available  concerning 
variations  in  geometry  and  properties.  The  probabilistic  solu¬ 
tion  relies  heavily  on  sensitivity  analysis  techniques,  and  for 
this  reason  is  applicable  to  models  which  are  large  and  complex. 
We  address  static,  steady-state  forced  vibration,  and  natural 
frequency  problems,  all  in  a  similar  fashion. 

Chapters  2  through  4  describe  the  analysis  methods  and 
solution  algorithms  used,  from  a  systems  point  of  view.  A  finite 
element  discretization  is  assumed,  but  the  development  is  other¬ 
wise  general.  We  begin  with  the  basic  solution  paths  in  Chapter 
2.  Chapter  3  develops  the  sensitivity  analysis  techniques  used 
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to  drive  the  probabilistic  calculations.  The  methods  described 
include  new  developments  in  geometric  shape  sensitivity  analysis 
which  also  have  potential  application  in  shape  optimization.  The 
probabilistic  analysis  is  then  outlined  in  Chapter  4. 

Chapters  5  and  6  deal  with  finite  element  technology.  In 
the  interest  of  streamlining  both  the  basic  solution  and  the 
sensitivity  calculations,  we  employ  a  Mindlin  plate  element  based 
upon  uniformly  reduced  numerical  integration.  While  substantial 
efforts  have  been  devoted  to  developing  improved  elements  of  this 

9 

type,  researchers  have  neglected  some  important  issues  which  are 
crucial  in  dynamic  problems;  this  is  the  main  topic  in  Chapter  5. 
Chapter  6  describes  the  methods  used  for  considering  layered 
components  using  conventional  plate/ shell  elements.  Although  the 
analysis  of  composite  blades  is  not  the  central  issue  in  this 
work,  it  is  likely  to  become  a  routine  requirement  in  the  future. 

Chapter  7  discusses  a  number  of  numerical  examples  which 
illustrate  various  aspects  of  the  present  study.  We  demonst *ate 
the  correctness  and  importance  of  the  element-level  techniques  of 
Chapters  5  and  6.  Several  sensitivity  analyses  show  the  capabil¬ 
ity  of  the  methods  described  here,  and  point  out  the  modeling 
techniques  which  are  preferred  for  geometric  sensitivities.  The 
probabilistic  solution  is  applied  to  analytical  examples  as  well 
as  comparisons  with  experimental  data. 

The  Appendices  contain  the  documentation  of  all  computer 
software  associated  with  the  work  described.  Computer  programs 
have  been  implemented  for  the  finite  element  solution  methods 
developed  herein,  and  for  data  communications  with  modeling  and 
graphics  software  such  as  PATRAN10  and  DISSPLA.11 


4 


CHAPTER  2 

ANALYSIS  PROCEDURES 


This  Chapter  reviews  the  basic  methods  of  solution  used  in 
the  deterministic  portion  of  the  finite  element  anlysis.  Primary 
emphasis  is  placed  upon  the  system-level  equations  resulting  from 
a  finite  element  discretization,  which  have  the  general  form: 

KU  +  MU  =  F  (1) 

Here  K  and  M  are  the  stiffness  and  mass  matrices,  which  normally 
are  large,  sparse,  and  symmetric;  U(t)  are  the  generalized  dis¬ 
placements  at  the  nodes  of  the  finite  element  model,  and  F(t)  are 
the  corresponding  generalized  forces.  The  details  of  the  finite 
element  approximation  are  described  separately  in  Chapters  5  and 
6.  For  more  information  on  the  basics  of  numerical  algorithms  as 
applied  to  finite  element  systems,  the  reader  is  encouraged  to 
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consult  the  standard  texts  on  finite  element  analysis. 


2.1  LINEAR  STATIC  SOLUTION 


For  quasi-static  loading  and  response,  the  applied  forces 
F(t)  are  constant,  and  the  accelerations  U  vanish.  The  system  of 
ordinary  differential  equations  (Equation  1)  describing  the  model 
then  becomes  the  algebraic  system: 

KU  =  F  (2) 

which  must  be  solved  for  the  (constant)  nodal  displacements  U. 

An  effective  solution  of  the  static  system  (2)  must  exploit 
the  symmetry  and  sparsity  of  matrix  K,  since  unique  nonzero  terms 
occupy  only  10-20  percent  of  the  matrix  in  most  problems.  The 
solution  technique  also  must  permit  economical  re-solution  of  the 
system  for  new  right-hand  vectors;  this  facility  is  useful  for 
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considering  multiple  loading  conditions  and  for  performing 
sensitivity  analysis. 

In  the  present  work  we  adopt  the  triple  factorization,  or 
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Gauss-Doolittle,  technique.  '  The  first  step,  which  involves 
only  the  coefficient  matrix,  is  the  symmetric  factorization: 

K  =  LDLT  ( 3 ) 


in  which  L  is  a  unit  lower  triangular  matrix: 


{  1.  i  =  j 

0,  i  <  j 


(4) 


and  D  is  diagonal.  It  is  straightforward  to  show  that  the  local 
bandwidth  of  L  never  exceeds  that  of  K.  Therefore,  the  strict 
lower  triangle  of  L  can  replace  that  of  K,  and  D  can  be  stored  on 
the  diagonal  of  K  to  minimize  storage  requirements.  Computations 
on  elements  outside  the  envelope  of  nonzero  coefficients  are  easy 
to  eliminate,  leading  to  an  efficient  solution  procedure. 


Once  the  factorization  is  complete,  the  equivalent  system 


LDLU  =  F 


(5) 


can  be  solved  in  three  steps: 


Lz  = 

F 

(forward  substitution) 

oy  = 

z 

(scaling) 

(6) 

ltu  = 

y 

(backward  substitution) 

for  each  right-hand  side  of  interest. 


The  equation-solving  routines  used  in  this  work  are  based 
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upon  the  solution  package  published  by  Felippa.  The  original 
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code  is  efficient  and  we 11 -documented,  and  has  been  tested 
extensively. 

2.2  NATURAL  FREQUENCY  SOLUTION 

In  the  natural  frequency  problem,  the  applied  forces  are  set 
to  zero  and  all  displacements  are  assumed  to  vary  sinusoidally  in 
time: 


U  =  Xsin(wt)  (7) 

The  resulting  discrete  equations  of  motion  become: 

KX  =  XMX  (8) 

2 

in  which  \  =  w  .  This  symmetric  generalized  eigenvalue  problem 
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is  solved  using  the  subspace  iteration  algorithm. 

Subspace  iteration  is  a  vector  iteration  method  in  which  a 
relatively  small  number  of  trial  vectors,  which  are  modified  in  a 
systematic  manner  to  span  the  least-dominant  p-dimensional  sub¬ 
space  of  K  and  M  (where  'p'  is  the  number  of  trial  vectors  used 
in  the  calculation) .  The  number  of  trial  vectors  is  selected 
automatically;  for  a  system  of  order  N,  for  which  n  eigenvalues 
are  to  be  computed,  we  take: 

p  =  min  [  N,  2n,  n+8  ]  (9) 

The  essential  steps  in  the  algorithm  are  as  follows: 

1.  Let  k=l,  and  define  starting  vectors  YQ. 

2.  solve  KXk+1  =  Yk  for  X 

I 

3.  Form  subspace  stiffness  =  x£+1KXk+1  =  *£+1*k. 
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r 


4.  Compute  Y^+1  =  MXk+1* 

T  T 

5.  Form  subspace  mass  matrix  »}c+1  =  X)C+iMX)c+i  =  Xk+lYk+l’ 

6.  Solve  the  subspace  eigenvalue  problem  (of  order  p) : 

*k+lqk+l  =  ®k+iqk+lAk+l 

for  the  diagonal  matrix  of  eigenvalues  A  and  the  eigen¬ 
vectors  q. 

7.  Arrange  the  eigenvalues  A  in  ascending  order;  normalize 
the  eigenvectors  q. 

8.  Form  new  trial  vectors  Y^+1  =  Xk+lqk+l* 

9.  Check  for  convergence;  for  each  eigenvalue,  convergence 

is  declared  whenever  4 — <  e  ,  where  6  is  a  toler- 

Ak+1 

ance  on  the  order  of  10~6. 

10.  If  not  converged,  let  k  «-  k+1  and  return  to  Step  (2)  . 


The  strong  points  of  the  algorithm  are  its  efficiency  for  large 
systems,  and  its  ability  to  maintain  a  respectable  convergence 
rate  for  systems  having  repeated  roots. 

For  unconstrained  systems,  the  stiffness  K  is  singular,  and 
the  solution  indicated  in  Step  (2)  of  the  algorithm  is  undefined. 
In  such  cases,  we  employ  an  eigenvalue  shift  which  renders  the 
coefficient  matrix  positive  definite.  In  place  of  the  original 
system,  we  solve; 


(K+sM) X  =  (X+s) MX 


(10) 
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for  the  shifted  eigenvalues  A+s  and  the  eigenvectors  X.  The 
shift  s  is  a  positive  number  which  must  be  sufficiently  large  to 
make  the  coefficient  matrix  (K+sM)  numerically  non-singular.  The 
eigenvectors  X  are  unchanged  from  those  of  the  original  system, 

and  the  natural  frequencies  may  be  recovered  using  w  =  J\ ,  after 
subtracting  the  shift  s  from  the  computed  eigenvalues. 

The  individual  mode  shapes  X  are  normalized  so  that  the 
magnitude  of  the  largest  displacement  component  is  equal  to  one. 
Stress  data  obtained  from  the  eigenvalue  solution  are  computed 
from  the  normalized  mode  shapes,  and  indicate  only  the  relative 
magnitudes  for  each  mode. 


2.3  STEADY-STATE  HARMONIC  SOLUTION 


In  steady-state  harmonic  analysis,  the  nodal  forces  vary 
sinusoidally  in  time,  so  that: 

F  =  F0sin(«t)  (11) 

where  u  is  a  known  forcing  frequency.  For  an  undamped  elastic 
system  the  steady-state  solution  U(t)  is  sinusoidal,  and  in  phase 
with  the  forcing  frequency, 

U  -  0Qsin (wt)  (12) 

The  dynamic  equations  of  motion  reduce  to: 

(K-w2M)U0  =  Fq  (13) 

For  a  given  frequency,  then,  the  problem  resembles  a  linear 
static  system  and  may  be  solved  directly  for  UQ.  The  usual 
stress  recovery  procedures,  based  upon  the  amplitudes  UQ,  result 
in  stress  amplitude  data,  since  a  =  a0sin(wt). 
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In  practice,  we  are  normally  interested  in  the  response  of 

the  system  throughout  a  specified  range  of  forcing  frequencies. 

The  solution  accepts  a  series  of  forcing  frequencies,  recomputing 
.  2 

the  harmonic  stiffness  K-u  M  at  each  frequency.  The  resulting 
displacement,  strain,  and  stress  amplitudes  at  selected  nodes  or 
elements  may  be  plotted  versus  forcing  frequency  to  characterize 
the  frequency  response  of  the  system. 
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CHAPTER  3 

SENSITIVITY  ANALYSIS 


This  chapter  describes  the  calculation  of  sensitivity  data 
by  direct  methods  for  isoparametric  plate  or  shell  elements. 
Sensitivity  parameters  of  interest  include  intrinsic  properties 
such  as  material  modulus  and  plate  thickness,  as  well  as  geometry 
variables  which  influence  the  size  and  shape  of  a  structure.  The 
sensitivity  calculation  therefore  must  consider  the  parametric 
mapping  within  an  element,  as  well  as  the  influence  of  geometric 
variables  on  the  orientation  of  an  element  in  space.  The  methods 
presented  specialize  directly  to  continuum  elements,  in  which  the 
coordinate  transformation  is  omitted,  or  to  simple  structural 
members  situated  arbitrarily  in  space. 

We  begin  with  the  development  of  the  general  relationships 
needed  for  performing  geometric  (or  shape)  sensitivity  analysis 
with  isoparametric  finite  elements.  The  additional  contribution 
to  geometric  sensitivity  caused  by  a  changing  local -to-global 
axis  transformation  is  considered  in  Section  3.2.  Finally,  the 
application  of  these  methods,  as  well  as  standard  techniques  for 
computing  property  sensitivities,  to  plate  and  shell  finite 
elements  is  discussed  in  Section  3.3. 


3.1  SENSITIVITY  FORMULAS  FOR  ISOPARAMETRIC  ELEMENTS 

This  section  presents  the  development  of  several  analytical 
relationships  needed  for  shape  sensitivity  calculations.  Methods 
for  computing  sensitivities  with  respect  to  intrinsic  properties 
of  an  element,  such  as  thickness,  density,  or  modulus,  are 

relatively  straightforward;  techniques  of  this  type  are  used 
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widely  in  structural  optimization.  When  control  parameters 

affect  the  nodal  positions  within  a  model,  however,  the  effect  of 

changing  a  given  parameter  is  much  more  complex.  Both  the  shape 

and  orientation  of  an  isoparametric  element  depend  upon  the  nodal 
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positions,  and  the  sensitivity  analysis  must  account  for  each 

effect  properly.  A  number  of  researchers  have  addressed  the 

23-27 

topic  of  geometric  sensitivity  analysis,  but  the  formulation 

of  efficient  computational  techniques  remains  an  important  area 
of  research. 

The  techniques  discussed  here  for  geometric  sensitivity 
analysis  are  oriented  toward  two-  and  three-dimensional  continua 
modeled  with  isoparametric  finite  elements.  The  notable  feature 
of  these  formulas  is  their  simplicity,  which  leads  to  quick  and 
systematic  computational  algorithms  for  most  standard  elements. 
While  our  use  of  these  methods  is  in  probabilistic  analysis,  the 
same  approach  is  suitable  for  use  in  structural  optimization. 

3.1.1  Shape  Functions  and  the  Jacobian  Determinant 

In  isoparametric  finite  elements,  element  stiffness  and  mass 
matrices  and  consistent  load  vectors  are  computed  by  numerical 
integration.  For  example,  the  element  stiffness  has  the  general 
form 


K  =  |n  bTDB  |j|  dne  (14) 

e 

in  which  B  contains  the  strain-displacement  relationship,  D  the 
elastic  constants,  and  |j|  is  the  Jacobian  determinant.  The  area 
or  volume  element  dfle  refers  to  the  unit  (or  biunit)  square  or 
cube  in  parametric  coordinates.  The  element  geometry  enters  this 
calculation  through  the  strain-displacement  matrix  B,  which 
consists  of  Cartesian  derivatives  of  the  element  shape  functions, 
and  the  determinant  J  J | .  In  a  two-dimensional  continuum  element, 
for  instance,  the  portion  of  the  strain-displacement  relation 
pertaining  to  node  I  of  an  element  is: 


B 


I 


3NI/3x1  0 

0  3N!/ax2 

3N]./ax2  d^1/dx1 


(15) 
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in  which  Nj  is  the  shape  function  for  node  I.  The  Jacobian 
determinant  is: 


lJi.l 

3x^ 

X  ^ 

xiK  d(a 

where  x...  is  the  coordinate  x.  at  nodal  point  K  of  an  element. 

The  coordinates  £  in  Equation  (16)  refer  to  parametric  direc¬ 
tions  within  an  element. 

First  consider  the  derivatives  of  the  elements  of  B  with 
respect  to  the  nodal  positions;  these  have  the  form  d (Nj  n)/dxk][. 
We  begin  with  the  identity  JJ-1=  I,  which  can  be  expressed  in 
indicial  form  as  follows: 

d£a 

3x^  JTq  =  ^a,mxm,0  *  *a,mNK,0xmK  =  5or/S  (17) 

Here  lower-case  Latin  indices  refer  to  the  Cartesian  coordinate 
directions,  upper-case  indices  to  the  nodes  of  an  element,  and 
Greek  indices  to  the  parametric  coordinate  directions  of  an 
element.  The  summation  convention  is  used  here  for  all  three 
types  of  indices;  a  comma  indicates  partial  differentiation  with 
respect  to  the  coordinate  following.  Note  also  the  interpolation 
used  for  the  spatial  coordinate  x^  within  an  element,  in  terms  of 
the  shape  functions  NR(£)  and  the  nodal  coordinate  values  xmR. 

Since  Equation  (17)  must  hold  for  any  values  of  the  nodal 
coordinates, 


^a,mNK,^xmK} 


axkl(6a0> 


0 


(18) 


Because  the  nodal  positions  are  independent  of  one  another,  we 
have  =  ^km^IK'  and  Ec3uati°n  (18)  becomes: 


/  — * (  C  \  — 

m,0  3xkIKa,m' 


"*a,kNI,0 


(19) 
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(20) 


or,  since  x^  ^ n  5mn' 


ax^^a.n1  ~^or,k*0,nNI,0  ~^a,kNI,n 


kl 


For  the  derivatives  of  NT  _  with  respect  to  the  nodal  coordin- 

j  j  n 

ates,  then,  we  obtain: 


axkI(NJ,n)  NJ,a  dxkI(^a,n)  ~NJ,o^a,kNI,n 


or 


axUT(NJ,n}  =  "NJ,kNI,n 


(21) 

(22) 


kl 


In  general,  if  the  nodal  coordinates  depend  upon  a  set  of  control 


parameters  P  ,  it  follows  that 


3x, 


dP  (NJ,n)  ~NJ,kNI,n  3Pmk 

m  n» 


(23) 


The  necessary  computations  to  determine  dB/dP^,  then,  involve 
only  the  original  shape  function  derivatives  and  known  data 
describing  the  dependence  of  the  nodal  coordinates  upon  the 
parameters  P 


m 


For  the  derivatives  of  j| ,  we  note  that 


J  =  X.  >  X 


ijfc 


(24) 


or  in  terms  of  the  nodal  coordinates: 


=  f ijkxiMxjNxkPNM,£ 1NN,C2NP,C 


(25) 


in  which  i  ...  is  t-he  permutation  tensor.  Note  that  each  nodal 
k 


coordinate  xk^  appears  linearly  in  Equation  (25);  it  follows  that 


and  evaluat- 


dljJ/dx^j  may  be  obtained  by  replacing  xk  Q  by  Nj 
ing  the  resulting  determinant.  The  expressions  so  obtained 
correspond  to  the  determinants  encountered  in  solving  the  system 


^a,kNI,k 


=  N 


for  Nj  ^  by  Cramer's  rule;  we  observe  directly  that 

SF  -  w  k 

£JXki  1,K 


(26) 


(27) 


Therefore,  the  derivatives  of  the  Jacobian  determinant  can  be 
computed  directly  using  only  the  shape  function  derivatives  and 
the  Jacobian  determinant  for  the  original  element.  When  the  node 
coordinates  in  turn  depend  upon  geometric  parameters  P  ,  we  have 


aPm 


dx. 


m 


(28) 


The  relationships  (22),  (23),  (27),  and  (28)  above  describe 

completely  the  dependence  of  the  element  matrices  upon  the  nodal 
positions,  and  provide  the  basis  for  many  important  sensitivity 
calculations.  The  next  two  subsections  illustrate  their  use  in 
some  common  cases  of  geometric  sensitivity  analysis. 


3.1.2  Stiffness  and  Stress  Sensitivities 

The  simplest  and  perhaps  most  common  sensitivity  calculation 
involves  the  determination  of  static  response  derivatives  with 
respect  to  control  parameters.  In  what  follows,  we  will  denote 
the  response  derivative  with  respect  to  a  typical  parameter  P  by 
a  prime;  that  is,  (  )'  =  d(  )/dP. 


For  the  displacement  sensitivities,  we  begin  with  the  static 
equilibrium  equations  Ku  =  F  for  the  complete  model,  and  note 
that 


K'u  +  Ku'  =  F' 


(29) 
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The  sensitivities  u'  therefore  may  be  found  using  the  already- 
factored  stiffness,  since 

Ku'  =  F'  -  K'u  (30) 

The  internal  force  sensitivity  K'u  in  Equation  (30)  is  best 
evaluated  directly  in  vector  form,  element  by  element,  and  then 
assembled  in  the  same  way  as  the  element  loading  vectors.  Since 
the  element  stresses  are 


a  =  DBu  (31) 

the  product  K'u  is  simply 

K'U  =  L  [(B')Tct|j|  +  BTDB'u|j|  +  BTa | J | ' ]  dfl  (32) 

e  e 

The  computation  of  K'u  is  possible  only  after  the  basic  solution 
is  complete,  but  requires  much  less  arithmetic  than  the  element 
matrices  themselves.  At  each  sampling  point,  it  is  necessary  to 
form  B  and  B',  compute  the  stresses  a  =  DBu  and  the  derivatives 
DB'u;  two  matrix-vector  products  then  complete  the  contribution 
to  the  integral,  since  the  last  two  terms  may  be  combined. 

Once  the  solution  for  u'  from  Equation  (30)  is  complete,  the 
stress  sensitivities  may  be  obtained  from 

O'  =  DB'u  +  DBu'  (33) 

for  each  element. 

3.1.3  Applied  Loads  Sensitivities 

The  load  sensitivity  F'  in  Equation  (30)  may  be  zero,  as  in 
the  case  of  point  forces,  or  may  depend  upon  the  model  geometry, 
as  for  pressure  loads  and  body  forces.  Consider  a  nonuniform 
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body  force,  whose  nodal  values  within  an  element  are  ;  that 
is,  the  components  of  the  force  vector  per  unit  volume  at  any 
point  are  Njf^j.  The  consistent  force  vector  is  then 

FkJ  -  J<)  fklNINJ  lJl  dne  <34> 

e 

Only  Jjj  is  affected  by  the  nodal  coordinates,  and  therefore 

fjLt  ‘  Jfl  fkININJ  lJl'  <35> 

e 

in  which  the  derivative  |j|'  may  be  computed  from  Equation  (28). 
Notice  that,  because  the  load  sensitivity  does  not  depend  upon 
the  displacement  solution,  and  F£j  may  be  computed  simul¬ 

taneously  with  the  original  consistent  loads  vector. 

Surface  load  sensitivities  are  more  complex,  since  geometry 
changes  may  affect  not  only  the  element  of  surface  area  but  the 
orientation  of  the  surface.  The  original  force  vector  involves 
the  surface  integral 


Fkl  -  i»  PNlnk  lJJ  d“e  <36> 

6 

Here  n)JJwldue  =  nkdAe  is  the  element  of  surface  area  in  physical 

coordinates,  and  o>e  refers  to  the  loaded  surface  in  parametric 

space.  As  for  the  body  forces  above,  the  pressure  can  be  inter¬ 
polated  from  nodal  pressure  values,  p({)  =  NK(£)pK.  If  R  denotes 

the  position  vector  of  a  point  on  the  surface  in  question,  then 

"  dAe  -  x  If;  a(idi2  <37> 

where  £  denote  the  parametric  coordinates  within  the  surface. 

In  terms  of  the  nodal  coordinates,  then, 
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Fkl  *  J’«e  P  e  i j kxiMx jN  "M.ijVfj  ni  d“e  <38> 

The  corresponding  derivatives  with  respect  to  nodal  positions  are 
given  by 

d  F 

NI  d"e  <38> 

Recall  that  the  derivatives  R  *  and  R  *  are  the  surface 

/?1  ><2 

tangents  needed  for  the  original  surface  pressure  calculation;  in 
terms  of  these  two  vectors,  formula  (39)  becomes; 


8  F, 


(40) 


PNitNJ,e1<VR,e2»-Nj,f2<VK,{;i»i  d“e 
or,  in  terms  of  a  design  parameter  P, 

f  5XnT 

Fkl  -  Jwe  PNI^J^1(inxR,e2)'NJ ((2(VR,f^lTd“e  <41> 


Again,  the  necessary  computations  depend  upon  quantities  which 
must  be  evaluated  to  form  the  original  load  vector,  and  can  be 
performed  at  the  same  time  as  the  consistent  loads  calculation. 


3 . 2  ORIENTATION  SENSITIVITY 


The  dimensionality  of  truss,  beam,  membrane  and  shell  finite 
elements  is  often  less  than  that  of  the  global  coordinate  system, 
and  element  calculations  must  be  performed  in  local  coordinates. 
Statistical  parameters  or  design  variables  which  affect  the  nodal 
coordinates  in  such  elements  control  both  the  element  dimensions 
and  orientation,  and  element  design  sensitivity  calculations  must 
account  for  both  effects.  The  influence  of  geometric  variables 
may  be  separated  into  two  distinct  contributions,  which  must  be 
applied  at  separate  stages  of  the  sensitivity  calculation. 
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In  components  built  up  from  bars,  beams,  and  panels,  or  in 
shell  structures,  one  additional  complication  arises.  Element 
stiffness  and  mass  properties  must  be  formulated  in  a  local 
coordinate  system,  and  then  transformed  to  a  common  coordinate 
system  for  assembly  and  solution.  Geometric  parameters  which 
affect  the  global  coordinates  X  at  a  node  may  influence  both  the 
element  coordinates  x  and  the  axis  transformation  relating  the 

two.  For  instance,  given  the  derivatives  for  a  single  para¬ 
meter  p,  and  the  global-to-local  transformation  x  =  AX, 


3x 

ap 


(42) 


The  effect  of  parameter  p  upon  A  cannot  be  neglected;  nor  can  the 
relationship  (42)  always  be  applied  directly  in  a  single  step. 

The  example  in  the  next  section  illustrates  both  of  these  points. 


In  subsequent  sections,  we  propose  a  simple  form  for  a 
general  coordinate  transformation,  for  which  derivatives  may  be 
computed  explicitly.  The  appropriate  calculations  are  outlined, 
and  the  method  is  applied  to  a  quadrilateral  element  in  three 
dimensions. 


3.2.1  Example  of  Orientation  Sensitivity 


The  planar  truss  problem  shown  in  Figure  2  demonstrates  the 
need  for  including  the  effect  of  geometric  parameters  on  the 
coordinate  transformation  for  an  element,  and  helps  to  clarify 
the  proper  methods  for  introducing  this  effect  into  the  sen¬ 
sitivity  calculation.  For  the  axial  force  member  in  the  Figure, 

dK 


we  wish  to  determine  the  derivative  =  JJ 


The  exact  result, 
for  K  referred  to  degrees  of  freedom  [UA, VA,U0, Vg) ,  is; 
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r  -2a& 

fl2  2 

-a 

2a$ 

a2-0 

£A 

(i2-<*2 

2a0 

.  2  -2 

a  -0 

-2a0 

L 

2a$ 

0 ,2->2 

-20tfl 

fi2-ct 

L  a2-fi2 

-2  a/5 

a2  2 

0  -a 

2a  0 

in  which  a=sin0,  fi-cos 8.  The  local  coordinate  transformation  is 

j"  x  j  =  cose  sind  ]  I"  X  ]  (44) 


cos#  sin0 
sin#  cos 8 


r x  ■ 

[  Y 


and  the  coordinates  at  end  '  B'  are  X0  =  Lcos0 ,  Yfi  =  Lsin# . 

d  A 

Since  8  does  not  affect  the  element  length,  neglecting  gj 

leads  to  x'  =  0,  and  hence  K'  =  0.  However,  it  is  easy  to  verify 
that  applying  equation  (42)  directly  yields  x'  =  y '  =0,  and  thus 
K'  =  0.  For  correct  results,  it  is  necessary  to  introduce  the 
two  contributions  in  equation  (42)  at  appropriate  stages  of  the 
element  calculation.  During  the  element  stiffness  computation  in 
local  coordinates,  only  the  overall  shape  and  dimensions  affect 
the  computed  results;  here  it  is  appropriate  to  introduce  the 

3  y 

"shape  effect",  ~  A jj,  holding  A  constant.  When  the  element 

matrices  are  transformed  to  global  axes,  the  "orientation  effect" 

=  |yX  is  significant.  If  the  local  stiffness  Kg  and  global 

stiffness  K  are  related  by 
9 


Kg  =  T  KgT 


the  appropriate  geometric  sensitivity  is,  in  general, 

*4  -  (T')TKjT  +  TTKp  +  TTKfT' 


in  which; 


d  K« 
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Here  3^  is  the  position  of  the  Kth  node  of  the  element  in  local 

coordinates.  The  range  of  summation  on  K  is  equal  to  the  number 
of  nodes  connected  to  the  element. 

The  observation  above  is  true  for  one-  or  two-dimensional 
elements  situated  in  three-dimensional  space  as  well.  In  prac¬ 
tice,  it  is  possible  to  avoid  much  of  the  computation  implied  in 
Equation  (46) ,  as  outlined  later  in  the  discussion. 

3.2.2  Basic  Coordinate  Transformation 


Below,  we  propose  a  form  of  the  local-to-global  coordinate 
transformation  which:  (a)  can  be  related  to  most  common  methods 
of  establishing  an  element  local  axis  system;  (b)  is  simple 
enough  that  analytical  expressions  for  its  derivatives  are  easily 
obtained;  and  (c)  requires  relatively  little  computation  to  form 
both  the  original  transformation  and  its  derivatives.  In  what 
follows,  we  assume  that  the  global  nodal  coordinates  depend  upon 
certain  geometric  parameters,  and  denote  a  typical  one  of  these 
by  "p" .  Furthermore,  the  effect  of  parameter  p  on  the  absolute 
location  of  the  element  centroid  is  neglected;  that  is,  we  let: 


?5s  _  ffs .  1  ?  !5i 

dp  dp  N  dp 


(48) 


in  which  N  is  the  number  of  nodes  per  element.  With  this  assump¬ 
tion,  the  origin  in  both  the  local  and  global  systems  may  be 
taken  to  coincide  at  all  times  without  loss  of  generality.  The 
effect  of  parameter  changes  on  absolute  position  must  be  account¬ 
ed  for  only  in  axisymmetric  elements,  for  which  the  coordinate 
transformation  is  often  unnecessary,  and  for  load  sensitivities 
which  depend  on  absolute  position,  such  as  centrifugal  forces. 


Consider  a  local  axis  system  defined  by  the  centroid  and  two 
additional  points  (Figure  3) .  The  positions  of  Points  1  and  2 
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relative  to  the  element  center  are  (X^Y^Z^  and  (X2,Y2,Z2), 

respectively.  The  local  x  axis  is  determined  by  the  centroid  and 
Point  1;  Point  2  provides  a  third  point  in  the  local  (x,y)  plane. 
In  terms  of  vectors  u  and  v,  shown  in  the  Figure,  the  unit 
vectors  defining  the  local  axes  are: 


u 


uxv 


el  ■  R!  e3  =  T^T  ®2  =  e3xel 


(49) 


Define  the  constants 


cyz  =  *1Z2 

"  Y2Z1 

czx  =  Z1X2 

‘  Z2X1 

Cxy  ~  X1V2 

-  X2Y1 

Dx  -  Zlczx 

"  YlCxy 

Dy  =  XlCxy 

ZlCyz 

Dz  =YlCyz 

-  XlCzx 

(50) 


(51) 


and  the  length  measures 


A  = 


=  yx| 

+  Y2  +  Z2 

-  J  c2 

+  C2„  +  C2  , 

yz 

zx  xy 

a  „  = 

o,a_ 

2 

1  3 

matrix  A, 

whose  rows  are 

simply: 

Xl/Qfl 

V°1  V°i 

Vtt2 

Va2  V°2 

CyZ/a3 

Czx/0r3  cxy/a3 

(52) 


(53) 
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3.2.3  Derivatives  of  Coordinate  Transformation 


Given  the  derivatives  X',  Y',  Z',  X^,  Y ' ,  Z',  it  is  a  simple 

matter  to  compute  the  derivatives  of  the  transformation  matrix  A. 
The  derivatives  of  the  constants  above  are  (for  example) : 


cyz  -  Yiz2  + 

Y1Z2  -  Y2Z1  - 

Y2Z1 

(54) 

Dx  =  ZlCzx  +  Z 

lCzx  -  YiCxy 

-  Yicxy 

(55) 

"  <X1X1 

+  Y1Y1  +  ZlZi 

)/«l 

(56) 

a '  =  (c  C' 

3  '  yz  yz 

+  C  c'  +  C  C'  ) /a_ 
zx  zx  xy  xy;/  3 

(57) 

*2  " 

aiQ3  +  °1°3 

(58) 

derivative 

of  A  is  then: 

Xlal-Xlal 

Y'a  -Y  a' 
ll  ll 

zial-zlai  1 

•? 

•i 

•l 

A'  = 

Di02-Dx°2 

Dv“2-Dv“2 

°;a2-w 

(59) 

q2 

(v  2 

a2 

°2 

a2 

a2 

c '  O'  -c  a ' 
vz  3  vz  3 

C'  Or  0 -C  Or ' 
zx  3  zx  3 

C'  a _ -C  a' 
xv  3  xv  3 

q2 

a  2 

/y  2 

L  °3 

°3 

a3  J 

It  is  easy  to  verify  that  the  original  transformation  A, 
computed  as  indicated  in  the  previous  section,  requires  10 
additions  or  subtractions,  28  multiplies  or  divides,  and  two 
square  roots.  If  the  derivatives  of  A  are  computed  at  the  same 
time,  an  additional  32  additions  or  subtractions,  61  multiplies 
or  divides,  and  no  additional  square  roots  are  required  per 
geometric  parameter. 
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3.2.4  Computational  Considerations 

In  practice,  computation  of  the  matrix  form  of  K'  is  usually 
unnecessary.  For  example,  sensitivities  for  static  analysis  may 
be  determined  from: 

Ku'  =  F'  -  K'u  (60) 

and  only  the  product  K'u,  formed  element-by-element  and  assembled 
as  a  vector,  is  required. 

Consistent  with  the  transformation  in  Equation  (45) ,  we 
assume  that  the  element  displacements  referred  to  local  coor¬ 
dinates  are  =  Tug.  Thus,  the  product  to  be  formed  is,  from 

equation  (46) : 

KgUg  =  (T')T(KgU£)  +  TT(Kg'u£  +  K?T'Ug)  (61) 

The  vector  in  the  first  term  represents  the  internal  element 

forces  in  local  coordinates,  which  may  be  computed  directly  from: 

=  Kluz  =  In  ®Ta  lJl  dae  (62) 

e 

The  vector  KgT'ug  appearing  in  the  last  term  can  be  obtained  in  a 

similar  fashion,  after  computing  the  stresses  corresponding  to  a 

fictitious  set  of  local  nodal  displacements  u^  =  T'ug.  As  noted 

in  Section  3.1.2,  an  efficient  means  of  calculating  the  remaining 
vector  KgUg  is  to  use  the  relationship 

Kg'Ug  =  Jn  [(B')Ta|j|  +  BTff'|j|  +  BTa  |  J  |  '  ]  dfle  (63) 

e 

Therefore,  the  calculation  of  sensitivities  related  to  the  local 
coordinate  transformation  requires  only  two  additional  internal 
force  evaluations,  and  transformation  of  the  resulting  vectors  to 
global  coordinates. 
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3.2.5  Example:  Line  Element  in  Space 


For  the  truss  member  considered  earlier,  and  the  geometric 
parameter  0,  Kg  -  0,  and  the  sensitivity  is  due  solely  to  the 

effect  of  0  on  the  coordinate  transformation.  Let  the  origin  of 
coordinates  correspond  to  Point  A,  and  Point  B  to  the  first  point 
defining  the  coordinate  transformation  (Point  1  above) .  Point  2 
may  correspond  to  any  point  in  the  plane  not  situated  on  the  line 
AB.  For  simplicity,  we  will  select  X2  =  0,  Y2  arbitrary.  Since 

X i  =  Lcos0 ,  =  Lsin0,  X'  =  -Lsin0,  Y'  =  Lcos0 ,  it  is  easy  to 

show,  from  Equation  (59),  that: 


A' 


-sin0  cos0  0 

-cos0  -sin0  0 

0  0  0 


(64) 


Using  Equation  (46)  with  Kg  -  0  yields  the  exact  result  shown  in 
Equation  (43) . 


3.2.6  Application  to  Plate  and  Shell  Elements 


Numerical  examples  for  a  two-dimensional  element  situated  in 
three  dimensions  are  somewhat  difficult  to  present  in  a  meaning¬ 
ful  form.  For  the  present,  we  simply  outline  the  application  of 
the  procedure  above  for  this  important  case.  Numerical  examples 
involving  orientation  sensitivity  are  presented  in  Chapter  7. 

Figure  4  shows  a  quadrilateral  element  in  three-space,  with 
a  common  choice  of  local  axes.  The  local  x  direction  is  oriented 
between  the  midpoints  of  edges  4-1  and  2-3;  a  vector  between 
midpoints  of  the  remaining  edges  completes  the  definition  of  the 
local  (x,y)  plane.  We  denote  the  coordinates  of  the  corners  by 
^XNi'  YNi' ZNi^  *  The  coordi-nates  and  x2  (for  example)  are  then: 
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X1  -  4  (-XNl+  XN2  +  XH3~  XK4> 

X2  =  4  *"XN1~  XN2+  XN3+  XN4* 

Coordinates  ,  y2,  Z2  are  defined  in  a  similar  way. 

derivatives  of  these  coordinates,  we  may  use: 


(65) 

For  the 


X<  =  I  <‘XN1+  XN2+  XN3-  XN4> 


N2 


N3 


X2  4  ^_XNl"  XN2+  XN3+  XN4 * 


(66) 


and  condition  (48)  is  satisfied  automatically.  At  this  point, 
Equations  (49)  through  (63)  apply  directly. 


3.3  SENSITIVITY  ANALYSIS  PROCEDURES 


The  preceding  Sections  present  the  mathematical  relation¬ 
ships  necessary  for  geometric  sensitivity  calculations  in  iso¬ 
parametric  finite  elements.  In  what  follows,  we  outline  the 
computational  procedures  used  in  sensitivity  solutions  for  a 
complete  finite  element  model.  The  sensitivity  parameters  of 
interest  include  intrinsic  variables  such  as  material  properties 
and  thicknesses,  as  well  as  geometric  control  variables  which 
govern  the  size  and  shape  of  the  model  by  controlling  the  nodal 
positions.  The  procedures  discussed  here  for  plates  and  shells 
may  be  specialized  to  other  isoparametric  elements  (where  the 
coordinate  transformation  is  omitted) ,  and  to  simpler  structural 
elements.  The  methods  described  are  efficient  and  accurate,  and 
relatively  simple  to  implement  for  most  standard  element  types. 

3.3.1  Element-Level  Calculations 


Consider  a  linear  static  problem  for  which  a  finite  element 
discretization  leads  to  the  algebraic  system  KU  =  F.  If  the 
stiffness  characteristics  of  the  system  or  the  applied  forces  are 
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dependent  upon  a  parameter  p,  then  the  dependence  of  the  nodal 
displacements  O  upon  p  may  be  obtained  by  solving: 


KU'  =  F'  -  K'U  (67) 

in  which  t  )'=d(  )/3p*  Notice  that  the  coefficient  matrix  in 
(67)  is  identical  to  that  of  the  original  problem,  so  that  the 
factors  of  K  may  be  reused  in  the  sensitivity  solution.  We  will 
focus  upon  the  calculation  of  the  product  K'U,  which  is  best 
performed  element-by-element,  and  then  assembled  for  the  complete 
system. 

While  it  is  possible  to  compute  K'  directly  for  an  element 
and  then  obtain  the  product  K'U,  this  approach  is  unnecessarily 
time-consuming.  We  prefer  to  form  K'U  directly  in  vector  form, 
which  reduces  both  the  number  of  arithmetic  operations  and  the 
computer  memory  required. 


Let  the  stiffness  matrix  for  an  element  be  given  by: 


K 


ATBTDBA  |j|  dQ 


e 


(68) 


in  which  A  is  a  transformation  from  local  to  global  coordinates, 
u=AU,  B  is  a  strain-displacement  matrix,  and  D  is  the  elasticity 
matrix.  The  region  is  the  domain  of  the  element  in  parametric 

coordinates.  The  transformation  matrix  A  may  vary  from  point  to 
point  for  curvilinear  elements,  but  is  constant  over  an  element 
in  most  simpler  elements. 


3. 3. 1.1  Intrinsic  Parameters 

The  product  K'U  is  simplest  to  obtain  when  parameter  p 
corresponds  to  an  intrinsic  property,  such  as  the  modulus  or 
thickness,  since  only  the  elasticity  matrix  is  affected.  Noting 
that  AU=u,  the  local  displacement  vector,  we  can  compute: 

0  =  D'Bu  (69) 
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and 


K'U 


T  T- 

A  B  ff 


|j| 


(70) 


The  computation  indicated  in  Equation  (69)  is  identical  in  form 
to  the  usual  process  of  stress  recovery,  so  that  the  calculation 
of  K'U  for  an  element  resembles  an  evaluation  of  the  internal 
nodal  forces.  In  fact,  the  element  internal  force  routines  may 
be  used  directly,  with  the  exception  of  calculating  D'. 


3. 3. 1.2  Geometric  Parameters 


When  the  parameter  of  interest  affects  the  nodal  posi¬ 
tions,  nonzero  derivatives  may  occur  for  the  element  of  area  | j| , 
the  strain  matrix  B,  and  for  the  coordinate  transformation  matrix 
A.  We  assume  that  the  derivatives  of  the  nodal  coordinates  are 
known,  and  represent  these  by  X£K=dX^K/dp,  in  which  i  ranges  from 

one  to  three,  and  K  from  one  to  the  number  of  nodes  per  element. 


The  calculation  of  B'  and  | j| '  depends  primarily  upon 

the  sensitivities  of  the  shape  function  derivatives,  d (N„  . )/dp 

1 

(Section  3.1.1).  The  transformation  sensitivity  A'  depends  only 
upon  the  global  nodal  coordinates  X^K  and  their  derivatives  XjK, 

as  discussed  in  Section  3.2.  From  the  relationships  developed  in 
Sections  3.1  and  3.2,  we  can  compute  the  product  K'U  from: 


K'U 


in  which: 


{  [(A')TBT+  AT (B')T]ff|j 


+  ATBT[ (o+o) I j|  +  a|j|']  }  dfte 

a  =  DBu 

a  =  DBu 

a  =  DB'u 
u  =  AU 

u  =  A'U 


(71) 

(72) 

(73) 

(74) 

(75) 

(76) 
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The  operations  indicated  in  Equations  (71-76)  are  analogous  to 
the  usual  displacement  transformation,  stress  calculation,  and 
internal  force  evaluation  steps  performed  in  a  linear  analysis. 


3. 3. 1.3  Mass  Matrix  Sensitivities 

Mass  sensitivity  calculations,  as  required  in  sen¬ 
sitivity  analysis  of  natural  frequencies,  are  simpler  in  form. 
Suppose  that  a  solution  has  been  performed  for  several  of  the 
dominant  modes  of  a  system: 

KU  -  (J2MU  =  0  (77) 

Differentiating  (77)  with  respect  to  the  parameter  of  interest 

28 

leads  to  the  frequency  sensitivity  expression: 

U?(K'-w?M')U. 

w?  =  - =r -  (i  not  summed)  (78) 

2wi0iM0i 

for  the  ith  mode  of  vibration.  Equation  (78)  remains  valid  when 
repeated  roots  are  present,  and  for  any  method  of  normalizing  the 
eigenvectors  IT.  The  denominator  is  a  scalar  multiple  of  the 

generalized  mass  for  mode  i,  which  we  choose  to  evaluate  at  the 

T 

system  level.  The  product  U^K'U^  may  be  computed  element  by 

element,  using  the  procedure  outlined  previously.  We  discuss  the 
evaluation  of  the  vector  M'U^  below. 


Letting  £,E  denote  a  particular  component  of  the  ele¬ 
ment  displacement  vector  in  local  and  global  axes,  respectively, 
we  write  the  contribution  to  the  mass  matrix  for  component  E  as: 


K-- 


/SATNNTA  |j|  dfle 


(79) 
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in  which  is  a  function  of  the  element  density  and  thickness. 

The  best  procedure  for  the  sensitivity  calculation  in  this  case 

T 

is  element  dependent.  However,  the  fact  that  N  As  =  £ (x) ,  the 
pointwise  value  of  £,  can  always  be  exploited.  Similarly,  the 
T  —  . 

product  H  A's,  resembles  a  point  displacement  value,  but  without 
the  same  physical  interpretation.  Again,  the  basic  sensitivity 
calculations  needed  are  limited  to  A/  and  | j| . 

3.3.2  System- Level  Solution 

This  section  outlines  typical  procedures  for  performing 

sensitivity  solutions,  assuming  that  element-level  routines  are 

available  for  evaluating  the  vectors  K'U  and  M'U,  and  the  scalar 

T  T  .  • 

products  U  K'U  and  U  M'U  as  required.  In  our  implementation  of 

these  methods,  we  perform  element  calculations  for  a  number  of 

load  cases  or  modes  and  for  a  number  of  sensitivity  parameters, 

all  in  parallel.  Sensitivity  parameters  may  include  the  material 

modulus  or  density,  element  thickness,  and  any  geometric  control 

parameter  defined  in  terms  of  derivatives  of  the  global  Cartesian 

coordinates  at  selected  nodes  with  respect  to  the  parameter. 

3.3. 2.1  Static  Response  Sensitivity 

In  static  analysis,  we  first  factor  the  original  stiff¬ 
ness  and  solve  for  the  nodal  displacements: 

K  =  LDLT  (80) 

LDLTU  =  F  (81) 

For  the  first  pass  of  sensitivity  calculation,  form  the  right 
hand  side  and  solve  for  displacement  sensitivities: 

Nel 

R  =  F'  -  l  (K'U)  (82) 

e  =1 
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ldltu' 


R 


(83) 


The  second  pass  of  element  sensitivity  calculations  yields  the 
element  stress  sensitivities: 

O'  =  (DBA) 'U  +  DBAU'  (84) 

Which  of  the  matrices  D,  B,  A  possesses  nonzero  derivatives  is  a 
function  of  parameter  type. 

Notice  that  if  a  particular  sensitivity  parameter  does 
not  affect  a  given  element  directly,  the  calculation  of  K'U  may 
be  skipped,  and  the  stress  sensitivity  reduces  to  a'  =  DBAU' .  In 
practice  it  is  convenient  to  maintain  a  list  of  switches  for  each 
element,  indicating  the  status  (active  or  inactive)  of  all  para¬ 
meters.  The  selection  of  parameters  such  as  modulus,  density, 
and  thickness  may  be  tied  to  material  or  property  set  numbers, 
making  it  easy  to  determine  whether  or  not  a  specific  element  is 
affected.  If  geometric  parameters  are  defined  in  terms  of  nodal 
coordinate  derivatives,  the  parameter  is  inactive  for  a  given 
element  only  if  all  derivatives  for  each  node  connected  to  the 
element  are  zero. 


3.3. 2. 2  Frequency  and  Mode  Shape  Sensitivity 


For  eigenvalue  problems,  we  first  solve  the  eigensystem 
and  compute  a  generalized  mass  for  each  mode: 


KUi  -  W^HO.^  =  0  (85) 

(i  not  summed) 

mi  =  U^MU.  (86) 

For  each  parameter  and  mode,  the  frequency  sensitivity  Equation 
(78)  may  be  summed  element  by  element: 


u 


i 


N 


2« 


el 

l 

e=l 


[  uTk'U 


0)?dTm'U.  ] 


(i  not  summed)  (87) 
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In  computing  sensitivities  of  the  eigenvectors,  we 
adopt  a  modal  representation,  as  suggested  in  Reference  28.  For 

the  i^*1  mode,  let  the  eigenvector  derivative  be: 

Uf  =  (88) 

in  which 

*  =  [  U1#  U2,  On  ]  (89) 

is  the  modal  matrix,  and  |l.  is  a  vector  of  modal  participation 

factors.  Introducing  (88)  into  the  derivative  of  the  original 

•  T  th 

eigenvalue  equation,  and  premultiplying  by  *  gives,  for  the  l 

mode: 

(k-w?m)^i  =  -tT(K'-w^M,)Ui  +  2wiu:#Tl«Ji  (i  not  summed)  (90) 

T  T 

Here  k  =  *Kt  and  m  =  #AMf  are  the  diagonal  generalized  stiffness 
and  mass  matrices.  Notice  that  only  the  ith  component  of  the 
product  t  MIT  ,  which  is  a  column  of  m,  is  nonzero.  The  element 
of  0^  corresponding  to  mode  'n'  is  therefore: 

-Un(K'_wiM')Di 

(0.)  =  - - -  (i*n;  u.tu  ;  i,n  not  summed)  (91) 

(k  -w.m  )  1  n 

'  nn  i  nn' 

Let  J  be  the  degree  of  freedom  which  attains  the  largest  value 
for  mode  i;  that  is: 

(U^j  =  sup  (U.)n  (92) 

2  8 

As  suggested  by  Rogers,  we  force  the  normalizing  basis  for  mode 
i  to  remain  constant  by  requiring  (U()T  =  0.  This  condition  is 

1  J 

sufficient  to  determine  the  remaining  element  of  0^: 
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(93) 


<*i>i  *  -  7ot)T  J,  l'i>»(0n>J 

i  J  n=l 

n*i 

in  which  N  is  the  number  of  modes  retained  for  the  sensitivity 

solution.  The  necessary  products  besides  the  system  generalized 

•  .  T  T 

stiffness  and  mass  consist  of  O  K'U,  and  U  M#U. ,  which  may  be 

n  i  n  i 

evaluated  on  an  element  basis. 

3. 3. 2. 3  Steady-State  Harmonic  Sensitivity 

The  steady-state  forced  response  solution  resembles  a 
linear  static  solution,  with  the  coefficient  matrix  K  replaced  by 
2 

K-w  M  for  a  given  forcing  frequency.  For  a  single  value  of  the 
forcing  frequency  u,  we  perform  both  the  basic  solution  and  all 
possible  sensitivity  solutions  together,  since  the  coefficient 
matrix  remains  constant.  The  numerical  procedure  is  precisely 
the  same  as  for  static  analysis,  with  obvious  changes  to  the 
coefficient  matrix  and  its  derivatives. 

The  interpretation  of  values  from  the  harmonic  response 
sensitivity  solution  is  different  from  static  or  free  vibration 
problems.  The  displacement  and  stress  solutions  now  represent 
amplitudes  of  these  quantities,  which  vary  sinusoidally  with  time 
at  the  forcing  frequency  (Section  2.3).  The  computed  sensitivity 
values  therefore  represent  derivatives  of  these  amplitudes  with 
respect  to  the  parameters  of  interest. 
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CHAPTER  4 

PROBABILISTIC  ANALYSIS 


This  chapter  outlines  a  general  approach  for  estimating  the 
variance  of  structural  response  variables,  given  the  mean  values 
and  variances  of  system  properties  which  are  probabilistic  in 
nature.  The  statistical  analysis  adopted  here  is  rudimentary,  to 
be  sure;  however,  the  treatment  is  consistent  with  the  level  of 
information  which  is  readily  available  to  the  engineer,  and  lends 
itself  to  the  analysis  of  relatively  large  and  complex  systems. 
The  sections  below  discuss  the  philosophy  of  the  approach,  the 
statistical  parameters  of  interest,  and  the  mathematical  rela¬ 
tionships  needed  for  computing  variances  of  response  quantities 
such  as  displacement,  stress,  and  natural  frequency. 

4 . 1  INTRODUCTION 


The  notion  of  a  probabilistic  analysis  encompasses  numerous 
possible  analytical  techniques.  Given  that  certain  properties  or 
dimensions  of  a  system  are  subject  to  uncertainty,  the  proper 
choice  of  analysis  method  depends  strongly  upon  the  difficulty  or 
cost  of  a  single  simulation,  and  upon  one's  knowledge  about  the 
statistical  parameters  of  interest.  We  note  some  of  the  possible 
approaches  below. 

Stochastic  analysis  involves  stating  the  differential  system 
of  interest  in  terms  of  stochastic  quantities,  and  solving  dir¬ 
ectly  for  the  response  in  statistical  terms.  This  approach  is  an 

.  .  .  29 

active  research  area  m  applied  mathematics  .  One-dimensional 

problems  still  represent  a  formidable  challenge  with  this  class 
30 

of  methods,  and  the  consideration  of  very  complex  systems  is 
not  feasible  at  this  time. 

Random  field  simulation  is  a  relatively  new  approach  devel- 
.  3 1 

oped  by  Liu  and  co-workers.  In  addition  to  discretizing  the 
deterministic  system  of  interest,  new  unknowns  are  introduced  in 
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a  finite  element  (or  other  numerical)  model  which  describe  higher 
statistical  moments  of  the  response.  This  augmented  problem  is 
solved  in  a  single  step  for  both  the  mean  values  (deterministic 
response)  and  the  additional  statistical  variables.  This  method 
is  capable  of  considering  detailed  autocorrelations  for  the 
statistical  variables,  but  is  best  suited  to  moderately  sized 
systems . 

Monte  Carlo  simulation  is  appropriate  if  the  statistical 
nature  of  the  (few)  independent  variables  is  well-known,  and  the 
cost  of  a  single  analysis  is  small.  Known  information  about  the 
statistical  parameters  is  used  to  generate  a  series  of  samples 
with  representative  values.  A  deterministic  analysis  is  per¬ 
formed  for  each  sample.  The  result  is  a  sample  of  the  response 

.  3  2 

from  which  statistical  data  can  be  derived  by  standard  methods. 

In  the  present  analysis  we  view  probabilistic  properties  of 
a  system  as  discrete  random  variables.  The  elastic  modulus  of  a 
turbine  blade,  for  instance,  might  vary  from  point  to  point  in  a 
different  fashion  for  every  blade  manufactured;  we  choose  to 
characterize  this  modulus  by  a  mean  value,  and  a  single  value  of 
the  variance.  The  information  needed  to  perform  a  meaningful 
probabilistic  analysis  with  this  approach  is  usually  available  or 
can  be  estimated  with  a  fair  degree  of  accuracy.  For  example,  if 
a  modulus  value  is  quoted  as  being  "EiAE’' ,  we  normally  interpret 
the  quantity  AE  as  representing  three  standard  deviations;  the 
range  EiAE  therefore  includes  approximately  99.7  percent  of  all 
samples. 

The  discrete  random  variable  approach  requires  a  similar 
level  of  information  about  all  statistical  variables.  Therefore, 
routine  quality  control  data  or  manufacturer's  tolerances  are  the 
only  additional  information  needed  beyond  that  used  to  construct 
a  deterministic  finite  element  model.  Together  with  the  rela¬ 
tively  low  cost  associated  with  solving  relatively  large  models, 
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this  simplicity  makes  the  present  method  attractive  for  routine 
analysis  work. 


It  should  be  noted  that  the  results  of  an  analysis  based  on 
the  discrete  random  variable  approach  are  not  related  in  a  simple 
way  to  results  from  the  alternative  methods  mentioned  previously. 
However,  the  basic  trends  predicted  by  either  method  will  agree; 
that  is,  if  the  dispersion  in  a  particular  variable  is  large,  and 
if  the  structural  response  varies  significantly  with  the  variable 
in  question,  then  we  expect  a  large  variance  in  the  response.  In 
some  cases,  it  is  possible  to  show  that  the  variances  predicted 
using  the  present  method  are  conservative  (overestimated) .  For 
example,  the  variance  in  natural  frequencies  predicted  when  a 
physical  property  is  assumed  to  vary  with  position  is  generally 
less  than  that  computed  when  the  property  is  constant  throughout 
the  model,  but  subject  to  the  same  variation  in  magnitude. 

4 . 2  STATISTICAL  PARAMETERS 

In  the  present  work,  we  consider  four  specific  types  of 
probabilistic  variables: 

□  elastic  modulus 

□  material  density 

□  thicknesses 

a  arbitrary  geometric  variables 

The  modulus  and  density  are  tied  to  a  specific  material  number  in 
the  finite  element  model.  Each  statistical  variable  is  defined 
by  specifying  a  property  set  number,  the  variable  type  (modulus 
or  density) ,  and  the  variance  expressed  in  consistent  units.  In 
a  similar  fashion,  a  thickness  variable  may  be  defined  for  any 
existing  property  set  in  the  model  by  specifying  the  number  of 
the  property  set  and  a  numerical  value  of  the  thickness  variance. 

Property  variables  (modulus,  density,  thickness)  represent 
the  simplest  cases  in  terms  of  finite  element  implementation, 
since  each  one  influences  the  stiffness  and  mass  characteristics 
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in  a  simple  way.  In  most  cases  the  stiffness  or  mass  matrix 
depends  linearly  upon  the  variable  in  question;  for  thickness 
variables,  bending  stiffnesses  vary  cubically;  however,  the 
necessary  computations  are  still  relatively  simple. 

"Arbitrary  geometric  variables"  are  implicitly  defined 
quantities  which  influence  the  overall  structural  geometry.  In 
the  context  of  a  finite  element  model,  these  geometric  variables 
influence  the  nodal  coordinates  for  all  or  part  of  the  model. 

The  geometric  variables  influence  the  stiffness  and  mass  charac¬ 
teristics  of  individual  elements,  but  in  a  more  complex  way  than 
for  property-based  variables. 


A  simple  example  of  a  geometric  statistical  variable  is 
useful  to  illustrate  the  nature  of  such  a  variable  and  the  tech¬ 
nique  used  to  define  it.  Figure  5  shows  a  segment  of  a  circular 
arc,  whose  radius  R  is  chosen  as  a  statistical  variable.  For  a 
node  located  on  the  arc  with  angular  coordinate  B  ,  the  Cartesian 
coordinates  of  the  node  corresponding  to  the  nominal  (mean)  value 
of  the  radius,  R^ ,  are; 

XH  =  Xc  +  Rficosff  f  Yn  =  Yc  +  Rnsin0  <94> 

In  specifying  the  nominal  coordinates  of  the  node,  the  value  of  R 
is  not  defined  explicitly.  We  define  the  statistical  variable  in 
terms  of  the  numerical  value  of  the  variance,  Var[R],  and  the 
effect  of  the  variable  on  the  existing  coordinates: 

§r  =  cos*  ;  §r  =  sin*  <95> 


For  a  response  variable  r(X,Y)  which  depends  upon  the  coordinates 

X  and  Y,  which  in  turn  vary  with  R,  it  is  possible  to  compute  the 
d  T 

derivative  without  knowing  the  nominal  value  of  R  explicitly; 


dj_  =  dJL  dX  ,  il  ix 

3R  (9X  dR  <9Y  dR 


(96) 
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Figure  5 


Circular  Arc  with  Variable  Radius. 


*4  1 


This  derivative  and  the  variance  of  R  are  sufficient  to  complete 
the  calculation  of  Var[T],  as  described  in  the  following  Section. 

4 . 3  VARIANCE  RELATIONSHIPS 


We  wish  to  formulate  a  relationship  between  the  mean  values 
and  variances  of  a  series  of  random  variables  (moduli,  densities, 
thicknesses,  and  geometric  parameters)  and  those  of  one  or  more 
structural  response  quantities.  Denote  the  random  variables  by 
p.,  and  a  typical  response  function  by  t (p1#p2 , . . . ,p  ) .  If  the 
function  r  is  linear  in  p^,  then: 

n 

r  =  l  a.p.  (97) 

i=l  1  1 

3  3 

then  the  expected  value  E(r)  is  simply 

n 

E[r]  »  l  a.E[p.]  (98) 

i*l 

and  the  variances  are  related  by: 


n  _  n-1  n 

Var[  t  ]  =  l  a^Var[  p.  ]  +  2  ^  £  a.a  .Cov[  p.  ,p.  ]  (99) 

i=l  i=l  j=i+l  J  J 

in  which  the  notation  Cov[a,b]  denotes  the  covariance.  When  the 
function  r  is  more  general,  linearization  of  r  about  the  mean 
values  [i  ^=e[  p^  ]  leads  to: 


Var[  r  ]  «  l  (fr-]zVar[Pi]  +  2  £  £  §T~fr-Cov[  p.,p .  ]  (loo) 

i=l  Fi  i=l  j  =  i+l  *1*3  x  J 

Note  that  the  derivatives  are  to  be  evaluated  at  the  point 

Pi^i '  i~l, 2 ,  •  •  • , n*  This  fact  is  exploited  in  constructing  an 
efficient  solution  procedure. 
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For  the  present,  we  will  consider  the  statistical  parameters 
of  interest  to  be  completely  independent,  so  that  Cov[p^,Pj]=0  in 
all  cases.  The  variance  of  any  response  quantity  r (p)  therefore 
is  given  by: 


vari')  -  ,l  [tk)2  Var[pi! 


(101) 


i=l 


Computation  of  the  response  variances  requires  only  the  variance 
of  each  statistical  parameter,  and  the  parameter  sensitivities 

3  T 

0^-.  Chapter  3  describes  these  calculations  in  detail. 


4.4  INTERPRETATION  OF  RESULTS 


The  results  of  a  probabilistic  analysis  include  a  basic 
solution  (expected  or  mean  values) ,  and  sensitivity  and  variance 
data  for  all  displacement  and  stress  variables.  Each  of  these 
components  of  the  solution  data  is  useful  in  its  own  way.  The 
basic  solution  identifies  the  type  of  response  which  occurs,  and 
is  used  to  identify  critical  areas  in  the  structure. 

Sensitivity  data  is  often  informative  for  design  purposes, 
because  it  indicates  which  variables  most  influence  the  response 
in  critical  regions.  It  is  important  to  recognize  that  the 
relative  magnitude  of  sensitivities  to  different  parameters  is 
not  necessarily  significant.  For  instance,  the  sensitivity  of  a 
displacement  or  natural  frequency  to  thickness  may  be  several 
orders  of  magnitude  larger  than  the  sensitivity  to  modulus,  only 
because  of  the  difference  in  magnitude  typical  of  these  two 
quantities.  In  many  cases,  the  product  ] /  where  a[  ] 

denotes  the  standard  deviation,  is  a  good  basis  for  comparing  the 
relative  importance  of  dissimilar  parameters. 

Statistical  data  generated  from  the  solution  are,  we  think, 
easiest  to  assimilate  when  presented  in  terms  of  a  single  scalar 
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quantity:  displacement  at  a  location  of  interest,  stress  at  a 
critical  point,  or  system  natural  frequency.  Variance  data  may 
be  of  interest  directly,  particularly  when  multiple  parameters 
are  involved.  Histograms  (bar-graphs)  are  often  a  useful  format 
for  presenting  and  assimilating  this  data,  because  they  reveal 
the  relative  influence  of  each  random  variable  on  the  uncertainty 
in  the  response. 


Another  concept  useful  in  interpreting  the  probabilistic 
analysis  results  is  that  of  a  percentile  value.  A  mean  value,  as 
computed  in  the  basic  finite  element  solution,  is  by  definition  a 
50th  percentile  value;  that  is,  we  expect  the  actual  value  to  be 
less  than  the  mean  in  50  percent  of  all  samples.  A  Q*"*1  percen¬ 
tile  value  Tq  is  such  that  the  true  value  of  the  variable  is  less 


than  or  equal  to  in  Q  percent  of  all  cases. 


Let  7  represent 


the  mean  value  of  the  variable  7,  and  7„  the  standard  deviation 
(the  square  root  of  the  variance) .  If  7  has  a  normal  (or  Gauss¬ 
ian)  probability  distribution,  the  interval  (7^-7^ , 7^+7^ )  repre¬ 
sents  about  68  percent  of  all  possible  values  of  7 .  The  interval 
(7^,7^+7a)  therefore  includes  approximately  34  percent  of  all 
possible  values  of  7;  this  means  that  the  value  of  7  will  be  less 

than  7  +7  84  percent  of  the  time.  That  is,  the  value  7..+r_  is 

£h  a  M  0 

the  84  percentile  value  of  7 .  The  percentile  value  also  can  be 

viewed  as  a  figure  of  reliability  or  confidence  level.  The 

changes  between  percentile  values  provide  a  direct  indication  of 

the  relative  uncertainty  in  a  particular  response  quantity. 


A  simple  example  is  useful  to  illustrate  the  interpretation 
of  a  probabilistic  solution  in  terms  of  percentile  values.  A 
thin  plate  supported  on  all  four  edges  has  modulus  E,  Poisson's 
ratio  v  ,  density  p ,  thickness  h,  and  side  length  a.  The  lowest 

4 

natural  frequency  of  the  plate  is  then:  4 

.2.  _ 

«  =  - § — ^  7E/p  (102) 

6  ( 1-1/ ^ )  a 
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The  sensitivity  of  this  frequency  to  the  plate  thickness  h, 
evaluated  at  the  nominal  thickness  value  hQ,  is: 

fh  ’  “'ho  <103> 


Accordingly,  if  thickness  is  the  only  statistical  variable,  we 
can  state  (from  Equation  101) : 


or 


var[y]  =  (M2  var[h] 

0 


(104) 


a  =  jr 
w  hQ  h 


(105) 


in  which  a are  the  standard  deviations  of  frequency  and 
thickness,  respectively.  Equations  (104)  and  (105)  apply  in  the 
neighborhood  of  hQ,  because  of  the  linearization  implied  in  (101) 
and  (104) . 


As  a  particular  case,  suppose  that  the  plate  is  made  from 

.  3  5 

2024T3  aluminum  flat  sheet,  mill  finish.  Manufacturer's  data 

for  actual  stock  thicknesses  indicate  that  acceptable  thickness 
variations  are  on  the  order  of  ±10  percent  for  very  thin  stock 
(less  than  0.030),  and  ±5  percent  for  thicker  sheet.  Interpret¬ 
ing  these  values  as  ±3ah,  we  take  ah=hQ/30  for  thin  sheet,  and 
ah=h0/60  for  thick  sheet.  From  Equation  (105) ,  the  standard 
deviations  of  the  natural  frequency  are  simply  ay=u/ 30  and 
ow=W/60  for  the  thin  and  thick  cases,  respectively. 

If  we  assume  a  particular  statistical  distribution  for  the 
statistical  variables,  we  can  interpret  the  result  in  terms  of 
percentile  values.  In  this  work  we  assume  a  normal  distribution 
for  all  variables.  For  a  normal  distribution  with  mean  n  and 
standard  deviation  a,  the  probability  of  a  value  less  than  or 
equal  to  'x'  is  given  by  the  distribution  function: 
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F(X) 


(106) 


aj2n 


r 

J-oo 


Vzit-|2 

21  a  J 


dv 


This  value  is  normally  written  in  terms  of  the  normalized  value 
z=(x-j u)/o;  in  effect,  the  "number  of  standard  deviations"  by 
which  x  is  separated  from  fi .  With  this  definition, 

Z  2 

F(x)  =  #(z)  =  y—  ln  e~U  /2du  (107) 

This  normalized  form  of  the  probability  distribution  function  is 
tabulated  in  most  statistics  texts  and  collections  of  mathemati¬ 
cal  tables.  The  meaning  of  i (z)  is  illustrated  in  geometric 
terms  in  Figure  6. 


For  a  0.99  confidence  level  (99th  percentile  value),  we  let 
$(z)=0.99,  and  from  the  statistical  tables  read  the  value  z=2.33. 
Recalling  the  definition  z =(x-n) /a ,  we  find  the  99th  percentile 
value: 


xgg  =  n  +  2.33a  (108) 

Values  of  the  normalized  variable  z  are  tabulated  for  selected 
percentile  levels  in  Table  1. 

For  the  plate  problem,  we  might  wish  to  determine  a  value  of 
natural  frequency  which  is  not  exceeded  by  most  plates.  To  bound 
Q  percent  of  all  cases  (the  percentile  value) ,  this  frequency 
is  Aw=w0+zaw,  where  z  is  the  value  of  (x-^)/a  corresponding  to  Q. 
In  other  words, 


1  +  z(Q) 


0) 


w. 


(109) 


Figure  7  shows  this  relationship  for  aw/wQ=  1/30  and  1/60;  the 
curves  are  labeled  "normal"  and  "high"  quality,  respectively. 
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Distribution  Function  $(z) 


7 


Table  1.  Number  of  Standard  Deviations  versus  Percentile  Level 


Percentile,  Q 

*(z) 

t 

z 

90. 

0.90 

1 . 282 

95. 

0.95 

1 . 6H5 

99. 

0.99 

2.326 

99.5 

0.995 

2.576 

99.9 

0.999 

3.090 

99.99 

0.9999 

3.719 

tx-u+zo,  with  z  corresponding  to  percentile  Q,  is  greater  than 
or  equal  to  the  sample  value  for  Q%  of  all  samples. 


Kiqurv  7.  Percentile  Values  of  Natural  Frequency 
for  a  Plate  with  Thickness  Variation. 


Note  that  the  99.9th  percentile  frequencies  are  about  12  percent 
and  6  percent  greater  than  the  nominal  natural  frequency,  due 
only  to  the  uncertainty  in  sheet  thickness. 
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CHAPTER  5 

FINITE  ELEMENT  APPROXIMATION 


This  Chapter  discusses  the  finite  element  approximations  to 
be  used  in  the  present  study.  The  choice  of  elements  is  dictated 
by  the  need  to  perform  accurate  solutions  for  both  thin  and  thick 
shells,  and  by  the  complexity  of  the  sensitivity  calculations.  A 
number  of  very  accurate  plate  and  shell  elements  exist  for  which 
the  sensitivity  computations  outlined  in  Chapter  3  become  hope¬ 
lessly  complex.  On  the  other  hand,  most  very  simple  elements, 
which  would  lend  themselves  to  compact  and  efficient  sensitivity 
computations,  do  not  possess  sufficient  accuracy  for  routine  use. 

The  finite  element  selected  for  use  in  this  work  is  a  four- 

node,  bilinear  displacement  element  based  upon  the  Mindlin  theory 
3  6 

of  plates.  Such  elements  exhibit  good  accuracy  for  both  thick 
and  thin  plates  when  reduced  (one-point)  numerical  integration  is 
used  to  evaluate  the  element  matrices.  However,  the  resulting 
element  is  rank-deficient,  and  must  be  "stabilized”  to  achieve 
reliable  behavior.  Methods  for  achieving  full  rank  of  the  stiff¬ 
ness  and  for  stabilizing  element  behavior  in  static  analysis  and 
in  explicit  dynamic  calculations  exist  and  are  quite  effective. 

To  date,  however,  very  little  attention  has  been  devoted  to  the 
proper  formulation  of  such  an  element  for  vibration  analysis  or 
in  implicit  transient  solutions. 

After  describing  the  origin  and  formulation  of  the  basic 
Mindlin  element,  this  Chapter  addresses  the  issue  of  controlling 
spurious  modes  of  response  in  dynamic  analysis.  Several  alterna¬ 
tives  for  the  element  mass  formulation  are  examined  in  detail. 

We  show  that  non-physical  dynamic  modes  exist  and  present  a 
potential  problem  with  most  mass  matrix  formulations,  and  that 
spurious  modes  other  than  the  familiar  hourglassing  motion  are 
possible.  A  combination  of  projection  methods  and  reduced  in¬ 
tegration  is  suggested  which  eliminates  these  deficiencies  and 
produces  accurate  numerical  results.  The  remaining  techniques 
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investigated  give  rise  to  anomalous  behavior  which  make  them 
unsuitable  for  general  use. 

5.1  BACKGROUND 

The  quadrilateral  Mindlin  plate  element  with  bilinear  dis¬ 
placement  and  rotation  fields,  based  on  single-point  quadrature, 

.  37 

was  introduced  by  Hughes,  Cohen,  and  Haroun  and  designated  ui. 

The  attractiveness  of  such  an  element  stems  from  its  simplicity, 

computational  efficiency,  and  high  accuracy  (since  the  single 

.  .  13 

quadrature  point  is  an  optimal  sampling  point  ) .  However,  the 

basic  Ul  element  is  rank-deficient,  since  bilinear  contributions 

to  the  displacement  field  are  not  captured  by  the  single-point 

integration.  Therefore,  the  assembled  stiffness  for  a  mesh  of  Ul 

elements  may  exhibit  singularities  when  properly  constrained,  or 

lead  to  the  prediction  of  spurious  oscillatory  displacements  with 

which  little  or  no  strain  energy  is  associated. 

Subsequent  development  of  the  bilinear  Mindlin  plate  element 
focused  largely  upon  the  stabilization  of  these  spurious  modes  of 
behavior.  In  the  context  of  explicit  dynamic  computations,  the 
concept  of  hourglass  stabilization,  as  discussed  by  Kosloff  and 
Frazier38  and  further  developed  by  Belytschko  and  co-workers39'40 
is  an  effective  means  of  controlling  this  behavior.  However,  the 
explicit  solution  provides  an  opportunity  for  individual  elements 
to  "react"  to  unstable  oscillatory  motions,  while  a  static  or 
implicit  dynamic  solution  does  not. 

41  4  2 

MacNeal  and  Hughes  and  Tezduyar  have  proposed  schemes 
for  stabilizing  the  bilinear  element  by  redefining  the  interpola¬ 
tion  of  the  transverse  shear  strain  field.  However,  these  tech¬ 
niques  require  a  four-point  quadrature,  and  the  simplicity  of  the 

A 

basic  element  is  lost.  Taylor^  and  Belytschko,  Liu,  and  co- 
44-47 

workers  have  pursued  the  idea  of  hourglass  mode  stabiliza¬ 

tion  for  static  analysis,  and  present  several  correction  methods 
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which  work  well  while  preserving  the  advantages  of  the  one-point- 
integrated  element.  Park,  Stanley,  and  Flaggs48  have  presented 
related  methods  of  stabilization,  obtained  as  a  by-product  of 
studies  on  element  behavior  with  increasing  mesh  refinement. 

Thus  far,  dynamic  calculations  based  upon  the  hourglass 
stabilization  techniques  appear  to  have  used  the  usual  consistent 
mass  for  the  quadrilateral  element,  and  very  little  information 
is  available  concerning  the  effect  of  the  stabilization  scheme  on 
dynamic  behavior.  Since  the  generalized  stiffnesses  used  in  the 
stabilization  scheme  must  be  small  to  avoid  locking,  one  suspects 
(correctly)  that  there  are  pitfalls  to  be  encountered  in  dynamic 
calculations.  Belytschko  and  Tsay  have  performed  eigenvalue 
solutions  for  plates  which  reveal  anomalous  behavior  in  some 
relatively  low-frequency  modes,  and  which  may  be  traced  directly 
to  the  effect  of  the  stabilization  operator. 

In  what  follows,  we  study  the  problem  of  formulating  the 
appropriate  mass  characteristics  for  the  bilinear  Mindlin  plate 
element  with  hourglass  stabilization.  The  purpose  of  this 
exercise  is  to  associate  the  proper  kinetic  energy  with  those 
motions  for  which  the  element  correctly  represents  strain  energy, 
and  to  eliminate  the  kinetic  energy  linked  to  the  spurious  modes. 
The  useful  frequency  range  should  be  free  of  vibration  modes 
controlled  by  the  stabilization  operator,  although  displacements 
associated  with  these  modes  must  be  permitted  to  occur  naturally 
as  a  part  of  lower-frequency  mode  shapes.  Modes  dominated  by  the 
stabilization  operator  should  be  relegated  to  the  higher  part  of 
the  frequency  spectrum,  which  is  already  dominated  by  the  finite 
element  discretization  rather  than  by  the  problem  physics. 

We  first  present  a  synopsis  of  the  bilinear  element  develop¬ 
ment,  and  introduce  some  useful  notation.  A  typical  stiffness 
stabilization  method  is  described.  We  then  examine  four  methods 
of  mass  matrix  formulation,  and  identify  their  most  important 
characteristics.  An  "optimum"  mass  formulation  is  suggested 
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which  associates  the  proper  kinetic  energy  with  basic  motions  of 
the  element  while  eliminating  spurious  dynamic  modes  caused  by 
the  stiffness  stabilization  scheme. 

5.2  BILINEAR  MINDLIN  PLATE  ELEMENT 

36 

The  kinematic  assumptions  of  Mindlin  plate  theory  relate 

the  displacements  (U,V,W)  at  a  generic  point  in  a  flat  plate  to 

displacements  (u,v,v)  and  rotations  (6  ,6  )  of  the  midsurface  by: 

x  y 

U(x,y,z)  =  u(x,y)  +  z*y(x,y) 

V(x,y,z)  =  v(x,y)  -  2*x(x,y)  (110) 

W(x,y, z)  =  w(x,y) 

in  which  z  is  the  direction  normal  to  the  midsurface.  The  state 
of  deformation  is  described  by  eight  generalized  strains, 

(t  -  [£x-V1,xy''VVltxy'Txz'1,yzI  (111) 

and  the  stress  state  by  the  corresponding  generalized  forces, 

at  =  [N  , N  , N  , M  ,M  ,M  , Q  ,Q  ]  (112) 

L  x'  y'  xy'  x'  y  xy'xxz  yzJ  v 

In  the  stress-strain  relationships  for  the  in-surface  strains  and 
curvatures,  plane  stress  assumptions  are  used.  The  transverse 

shear  quantities  are  related  by  Qqz  =  kGt702,  in  which  k  is  a 

•  49  •  •  • 

shear  correction  factor  ;  here  we  consider  only  isotropic  plates 

5 

and  employ  the  value  k-~  throughout. 

O 

For  the  bilinear  finite  element  (Figure  8) ,  we  use  the  shape 
functions: 

n  =  J  (  s  +  ££  +  nn  +  htn  )  (ii3) 

in  which: 

st  =  [  1,  1,  1,  1  ]  (114) 

=  [-1,  1,  1,-1  ]  (115) 
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n t  =  [-1.-1,  l,  i  ]  (lie) 

=  [  1,-1,  1,-1  ]  (117) 


47 

Following  Liu  and  Belytschko  ,  we  also  define  the  following 
useful  quantities: 


_1 

2A 

t*24'*31'y42'»13> 

(118) 

1 

2A 

'X42'*13'X24'X3ll 

(119) 

in  which  x..=x.-x.  and  y..=y.-y. .  Note  that  the  element  area  is 

-  Zj  1  j  lj  Z  J 

A=^(x31y42+x24y31) .  With  some  algebra,  it  is  possible  to  verify 
that: 


-  dN 
bl  = 


t=T)  =  0 


•  5N 
b2  =  3y 


e=r?=o 


(120) 


Finally,  we  will  make  use  of  the  element  corner  coordinates  in 
the  form  of  the  four-dimensional  column  vectors: 


x  =  [x1,x2,x3,x4] 


y  =  [y1,y2,y3,y4] 


(121) 

(122) 


It  is  useful  to  note  the  following  relationships  which  exist 

.  .  45 

among  the  vector  quantities  defined  above: 

sfch  =  stb1  =  sfcb2  =  htb1  =  hfcb2  =  bjy  =  b^x  =  o  (123) 


.  t  .  t 

bix  =  b2y  =  1 


(124) 


Equations  (123)  and  (124)  are  particularly  useful  in  identifying 
the  linearly  independent  modes  of  behavior  for  the  element.  For 
example,  u=s,  where  u  are  nodal  displacements,  defines  a  uniform 
(rigid-body)  motion,  while  u=h  defines  one  possible  hourglass 
deformation  pattern  in  a  rectangular  element  (Figure  9) . 
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Figure  9.  Hourglass  Displacement 
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With  the  definitions  above,  the  strain-displacement  operator 
evaluated  at  the  element  center  can  be  written  as: 


Writing  the  element  displacement  vector  as: 

dfc  =  [  u,  v,  v,  $x,  ]  (126) 

then  e=Bd  are  the  element  generalized  strains  sampled  at  the 
centroid  of  the  element.  With  the  stress-strain  relation  ff=De , 
the  element  stiffness  obtained  through  one-point  guadrature  is: 

K  =  I  BtDB  dA  =  BtDBA  (127) 

J  A 

5.3  STIFFNESS  MATRIX  STABILIZATION 

We  will  employ  a  stabilization  technique  for  the  stiffness 
matrix  based  upon  the  generation  of  hourglass  suppression  forces, 
as  suggested  in  References  44-47.  The  plate  element  of  the 
previous  section  contains  twenty  degrees  of  freedom;  we  first 
segregate  the  possible  modes  of  behavior  for  the  element  into  the 
eight  uniform-strain  modes  captured  by  the  single-point  quadra¬ 
ture  rule,  the  six  proper  rigid  body  motions,  and  the  six  remain¬ 
ing  modes  which  must  be  stabilized.  The  rigid  body  motions  are: 
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u  =  s 

(x  translation) 

V  =  s 

(y  translation) 

V  =  s 

(z  translation) 

(0 

II 

X 

>1 

II 

> 

(x  rotation) 

w  =  -X  9y  =  s 

(y  rotation) 

u  =  -y  v  =  x 

(z  rotation) 

The  six  spurious  inodes  consist  of  five  hourglass  deformation 

patterns  u=h,  v=h,  w=h,  9  =h,  0  = h,  and  the  twisting  mode: 

x  y 

W  =  J(sfcy)x  +  ^(stx)y;  $x  =  x;  0y  =  y  (128) 

The  stabilization  scheme  must  inhibit  the  hourglassing  modes  to 

produce  a  stable  element.  The  spurious  twisting  mode  exists  for 

a  single  element,  but  cannot  occur  in  a  mesh  of  two  or  more 

50 

elements,  as  shown  by  Hughes. 

4  5 

Belytschko  and  Tsay  define  generalized  hourglass  strains 

associated  with  each  component  of  displacement  and  rotation  by 
t  t 

7  and  7  $a ,  in  which: 

7  =  h  -  (hfcx) bx  -  (h^bj  (129) 


The  last  two  terms  of  Equation  (129)  are  important  for  irregular 
elements,  if  the  hourglass  strains  are  to  vanish  in  the  presence 
of  rigid  body  motion  and  uniform  strain.  We  will  use  a  defini¬ 
tion  which  is  only  slightly  different,  letting: 


=  2 


X3iy31*  X24y42 


7  = 


d  2N 


dxdy 


£=r)=0 


(130) 


In  the  stabilization  method  used  for  the  present  study,  we 

define  a  generalized  hourglass  strain  associated  with  each  of  the 

(h)  — t 

displacement  components  (e.g.,  f ^  -  7  u) /  and  associate  with 

this  strain  a  generalized  stiffness  determined  through  numerical 
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experiments.  The  generalized  hourglass  stiffnesses  for  the 
individual  displacement  components  are: 


E(h)  .  E(h) 

u  v 


0.10  Et 

(  1  +  A  > 


E(h>  ,  E(h,  ,  E(h) 

x  y 


0.10  Et* 


(  1  + 


A  > 


(131) 


(132) 


The  factor  of  (1+j^)  in  the  hourglass  stiffnesses  is  motivated  by 
locking  problems  observed  in  elements  with  extremely  small  dimen¬ 
sions.  For  small  element  areas,  the  correction  terms  contain  an 
additional  factor  proportional  to  the  area,  which  reduces  the 
artificial  strain  energy;  when  the  element  area  is  significant, 
the  1/A  contribution  becomes  small.  While  the  above  scheme  is 
not  optimal  for  some  problems,  we  have  obtained  good  behavior  for 
aspect  ratios  0 . 0001< (L/t) <0 . 1  over  a  range  of  six  orders  of 
magnitude  in  the  planform  dimension  L. 


5.4  EFFECT  OF  STABILIZATION  IN  DYNAMICS 

The  stabilization  scheme  outlined  above  is  effective  for 
static  problems,  in  which  the  anti-hourglass  stiffnesses  may  be 
adjusted  freely  to  produce  good  element  behavior.  In  dynamics, 
however,  low-energy  deformation  patterns  consisting  mainly  of 
hourglassing  motions  represent  likely  modes  of  low-frequency 
oscillation;  the  resulting  non-physical  solutions  may  contaminate 
that  portion  of  the  vibration  spectrum  which  frequently  is  of 
greatest  interest.  It  is  useful  to  view  this  problem  in  terms  of 
generalized  stiffness  and  mass  quantities,  as  follows. 

Given  a  vibration  mode  shape  d=t ,  the  generalized  stiffness 
and  mass  associated  with  the  mode  are  the  projections: 

k  =  *tK*  m  =  ♦tM*  (133) 
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The  corresponding  frequency  of  vibration  is  w=7k/m.  Since  the 
hourglassing  inodes  identified  in  the  preceding  section  are 
orthogonal  to  both  the  rigid  body  motions  and  the  constant  strain 
states,  the  generalized  stiffness  associated  with  hourglass  modes 
depends  solely  upon  the  hourglass  stiffnesses  E^,  which  are 
made  small  deliberately  to  avoid  locking  problems.  It  is  easy  to 
show  that  at  the  same  time,  the  kinetic  energy  associated  with 
the  unstable  modes  is  similar  in  magnitude  to  that  of  the  rigid 
body  and  uniform  strain  motions.  Consequently,  artificial 
vibration  modes  whose  frequency  is  governed  exclusively  by  the 
generalized  hourglass  stiffnesses  appear  low  in  the  element 
spectrum,  intermixed  with  the  lower-frequency  vibration  modes 
which  are  commonly  of  interest. 

This  difficulty  can  be  corrected  by  reducing  or  eliminating 
the  kinetic  energy  associated  with  unstable  modes  of  the  element, 
while  leaving  the  energy  associated  with  rigid  body  motions  and 
uniform  strain  states  unchanged.  In  other  words,  we  wish  to 
minimize  the  projection  of  the  mass  matrix  on  the  unstable  modes 
of  the  element,  so  that  spurious  modes  of  vibration  occur  only 
beyond  the  useful  frequency  range  of  the  model.  At  the  same 
time,  motions  which  are  legitimate  but  which  contain  hourglassing 
components  (such  as  torsional  or  inplane  bending  modes)  may  occur 
with  only  the  minimal  constraint  imposed  by  the  antihourglassing 
mechanism. 

5.5  MASS  MATRIX  FORMULATION 

In  this  Section,  we  examine  several  possible  constructions 
of  the  bilinear  element  mass  matrix.  The  point  of  this  exercise 
is  to  identify  a  mass  formulation  which  is  reliable  when  used  in 
conjunction  with  the  hourglass  stabilization  technique  described 
previously.  We  show  that  among  the  obvious  choices  for  the 
element  mass  formulation,  there  is  only  one  method  which  elimin¬ 
ates  the  possibility  of  unstable  solutions  for  dynamic  problems. 
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5.5.1  Fully  Integrated,  Consistent  Mass 


For  the  mass  properties  of  the  bilinear  element,  we  define: 


and 


ft/2  2 

( R- , R~ , R3 ) 

p(l,z,z  )  dz 

(134) 

1  ^  J 

J  -t/2 

H  = 

|  HNt  dA 

(135) 

A 


The  kinetic  energy 

T  =  2  J  [R1(u2+v2+w2)  +  2R2(u$y-v$x)  +  R3(*x+*y>]  dA  (136) 

A 

then  leads  to  the  consistent  mass  matrix: 


-  R^H 

0 

0 

0 

R2H  1 

0 

RjH 

0 

-R2H 

0 

0 

0 

RXH 

0 

0 

(137) 

0 

-R-H 

6 

0 

R-H 

6 

0 

L  r2« 

0 

R3H  . 

With  the  assumption  of  a  constant  Jacobian  determinant,  matrix  H 
may  be  evaluated  directly  in  closed  form,  giving: 


H  = 


2 

4 

2 

1 


1 

2 

4 

2 


2  1 
1 
2 

4  J 


(138) 


5.5.2  Lobatto  Integrated,  Consistent  Mass 

A  lumped  mass  matrix  for  the  bilinear  element  can  be  formed 
using  Lobatto  integration16,  with  the  quadrature  points  placed  at 
the  four  nodes  of  the  element.  The  resulting  matrix  is  diagonal 
provided  R2=0,  which  is  the  case  for  isotropic  plates  or  midplane 
symmetric  laminates.  The  translational  and  rotational  inertias 
associated  with  each  node  are  ^R^  and  4R3A,  respectively. 
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5.5.3  Consistent  Mass  via  a  Projection  Method 


One  obvious  problem  exists  with  the  consistent  mass  matrix 
evaluated  with  2x2  quadrature,  in  that  the  kinetic  energy  of  the 
hourglass  modes  is  relatively  large  compared  with  the  hourglass 
stiffnesses  used  for  stabilization.  Anticipating  this,  we  adopt 
a  simple  projection  method,  designed  to  eliminate  the  kinetic 
energy  associated  with  the  hourglassing  modes,  and  formulate  the 
mass  matrix  as  follows. 

Consider  the  kinetic  energy  associated  with  the  inplane 
displacement  u: 


u 


-  M 


Riu 


dA 


(139) 


We  wish  to  eliminate  the  kinetic  energy  associated  with  that  part 
of  the  velocity  field  ignored  by  the  single-point  integration  of 
the  element  stiffness.  To  this  end,  we  expand  the  velocity  field 
u(x,y)  about  x=y=0  (which  we  assume  coincides  with  the  element 
center) : 


U(x,y)  =  u  (0 , 0)  +  xu  v(0,0)  4-  yu  (0,0)  +  ...  (140) 

or 

u(x,y)  =  ^u  +  x(b^u)  +  y  (b^u)  +  £77  (^u)  (141) 

We  now  form  a  modified  kinetic  energy  based  on  the  purely  linear 
part  of  the  velocity  field, 

U£ (x,y)  =  Js1^  +  x(b^u)  +  yfb^u)  (142) 

giving: 

=  2  R1(4s+xb1+yfa2) (4s+xbi+yb2) dAl  «  (143> 

A 
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As  with  the  stiffness  computation,  we  assume  a  constant  Jacobian 
determinant  over  the  element.  Note  that,  using  (113), 

f  H  dA  =  f1  f1  j(s  +  ££  +  tit)  +  hfrj)  |  j|  d£dr?  =  |j|s  (144) 
J  A  J  -l 

If  the  center  of  the  element  is  x=y=0,  then: 

x  dA  =  f  dA  =  |  j|  (stx)  =  0  (145) 

A  J  A 

and  terms  linear  in  x  or  y  do  not  survive  the  integration.  Using 
(135),  we  define: 

IA  x*  dA  =  1't,,x  *  cxx 

lA  y2  dA  -  yt"y  -  cyy 

xy  dA  =  xtHy  =  cxy 

and  the  modified  kinetic  energy  becomes: 

T?  *  \  R1UtHu  (149) 

in  which: 

s  -  ^f-sst  +  oxxblbJ  +  cyyb2b‘  +  cxy(blb2+b2bl)  (150) 

The  evaluation  of  H  requires  no  numerical  integration,  and  the 
necessary  vector  products  are  identical  with  existing  terms  in 
the  element  stiffness  matrix. 

Performing  a  similar  linearization  of  all  displacement  and 
rotation  components,  we  find  that  the  projected  element  mass 
matrix  is  identical  in  form  to  equation  (137),  with  H  replaced  by 
H.  Using  Equations  (123)  and  (124),  it  is  evident  that  the 
kinetic  energy  (and  therefore  the  generalized  mass)  associated 


(146) 

(147) 

(148) 
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with  all  five  pure  hourglassing  modes  is  identically  zero?  there¬ 
fore,  inodes  consisting  primarily  of  hourglassing  motions  should 
not  appear  as  spurious  low-energy  vibration  modes. 

5.5.4  Consistent  Mass  by  Reduced  Integration 


With  a  single  point  quadrature,  the  integral  in  Equation 
(133)  is  sampled  only  at  the  centroid  of  the  element,  where  N=^s. 
The  resulting  mass  matrix  is  then  identical  to  that  of  Equation 
(137),  with  H  replaced  by: 


H 


(1) 


1  1  1 
1  1 

1  1 

1  1  J 


(151) 


Note  that  the  mass  matrix  so  obtained  should  be  equally  effective 
to  the  projection  method  in  eliminating  the  kinetic  energy  of  the 
hourglass  modes,  which  are  bilinear  in  x  and  y. 


5.5.5  Comparison  of  Mass  Matrix  Formulations 


Certain  properties  of  the  four  mass  matrix  formulations 
described  above  are  readily  apparent,  and  numerical  experiments 
(see  the  next  section)  reveal  additional  characteristics  which 
are  less  obvious.  We  discuss  the  most  important  of  these  below. 
In  all  cases  we  assume  a  stiffness  matrix  formed  using  single 
point  integration,  and  stabilized  using  hourglass  control. 

With  a  fully-integrated,  consistent  mass  (Equation  137) ,  the 
occurrence  of  spurious,  low-frequency  hourglass  modal  patterns  is 
expected.  The  kinetic  energies  associated  with  hourglassing  are 
similar  in  magnitude  to  that  of  other,  legitimate  vibration 
modes,  while  the  stiffnesses  are  typically  an  order  of  magnitude 
less.  The  generalized  stif fness-to-mass  ratios  for  hourglassing 
patterns  are  therefore  quite  low,  and  spurious  modes  will  appear 
low  in  the  vibration  spectrum. 
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Similar  problems  with  the  lumped  mass  formulation  are  to  be 
expected,  since  a  diagonal  mass  matrix  leads  to  kinetic  energy 
for  all  possible  motions.  Likewise,  any  positive  definite  con¬ 
sistent  mass  formulation  is  destined  to  predict  non-physical 
dynamic  motions. 

The  mass  matrix  obtained  by  projection  onto  a  linearized 
velocity  field  is  positive  semi-definite,  since  zero  kinetic 
energy  is  associated  with  all  pure  hourglassing  patterns.  This 
method  is  satisfactory  for  plate  bending  alone,  since  the 
singular  modes  of  the  mass  matrix  coincide  precisely  with  those 
of  the  stiffness.  Such  is  not  the  case  for  inplane  behavior, 
since  spurious  low-energy  modes  which  do  not  correspond  to  pure 
hourglass  patterns  may  still  occur. 


The  troublesome  vibration  mode  of  Figure  10  involves  hour- 
glassing,  combined  with  a  uniform  rotation  about  the  element 
centroid.  In  a  rectangular  element  of  dimension  (2a, 2b),  for 
example,  the  inplane  motion  can  be  described  by: 


-  o  - 

U  =  2  (*I+h)  = 

-  o  - 

v  =  -  5*(C+h)  =  -  § 

■  0  - 

l  (152) 

-  -1  - 

the  rectangular  element 

,  • 

■  • 

in  question, 

equations 

( 118 ) - ( 119 ) 

H  1  C  . 

1  _ 

(153) 

bl  "  4a  ?  ' 

b2  = 

4b  1 

Since  vectors  f  and  tj  are  orthogonal,  the  resulting  centroidal 
strains  vanish,  though  the  motion  is  not  a  rigid-body  rotation 
The  projected  mass  neglects  the  hourglass  velocity  components, 
and  predicts  a  kinetic  energy  based  on  the  nodal  velocities: 


u. 


(34) 
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corresponding  to  a  uniform  rigid  body  rotation  about  the  element 
center.  It  is  worth  noting  that  the  same  non-physical  velocity 
field  is  possible  at  low  frequency  with  the  exact  consistent  mass 
matrix,  though  the  kinetic  energy  is  higher  due  to  the  presence 
of  hourglass ing. 

Figure  11  shows  the  deformation  pattern  developed  within  a 
uniform  mesh  of  rectangular  elements  for  a  mode  of  this  type. 

Dark  lines  indicate  element  edges  which  experience  only  pure 
ex;ension,  compression,  or  rigid-body  translation. 

Pure  hourglass  motions  and  the  staircase  patterns  of  Figure 
4  are  both  made  possible  by  the  bilinear  element  displacement 
field.  The  hourglass  field  is  a  bilinear  displacement  field 
which  vanishes  at  the  element  center,  causing  zero  strain,  while 
the  staircase  mode  consists  of  a  constant  rotation  field  about 
individual  element  centers,  combined  with  hourglassing  motion. 

Correction  of  the  spurious  inplane  mode  problem  requires 
that  only  true  rigid-body  rotations  lead  to  a  nonzero  kinetic 
energy.  The  use  of  a  single-point  quadrature  achieves  this 
property,  since  kinetic  energy  results  only  for  mean  rotations 
about  a  point  other  than  the  element  centroid.  The  only 
unfortunate  consequence  of  this  choice  is  that  the  kinetic  energy 
of  an  element  whose  centroid  coincides  with  an  axis  of  inplane 
rigid  body  rotation  will  be  missed. 

Deformation  patterns  analogous  to  the  inplane  staircase  mode 
do  not  appear  to  exist  for  out-of-plane  vibration.  The  projected 
mass  matrix  therefore  may  be  used  with  confidence  for  flat  plate 
bending. 

It  remains  for  us  to  compare  the  single-point  integrated 
mass  matrix  with  the  projected  mass  (both  of  which  are  immune  to 
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hourglass  inodes)  in  the  bending  problems  for  which  both  are 
potentially  applicable.  Consider  the  out-of-plane  rigid-body 
modes  described  earlier,  and  the  uniform  strain  states: 

tfy  =  cx  (x  curvature) 

0x  =  cy  (y  curvature) 

*x  =  clx  %  =  c2y  (twist) 

w  =  c1x  Oy  =  c2s  (x-z  shear) 

w  =  c1y  =  c2s  (y-z  shear) 

One  useful  comparison  is  based  on  the  kinetic  energy  associated 
with  each  of  these  motions  using  the  consistent,  projected,  and 
reduced-integrated  mass  matrices.  We  find  that  the  projected 
mass  matrix  yields  the  proper  energy  for  all  eight  elementary 
states,  and  zero  for  the  hourglass  modes.  The  mass  obtained  by 
single-point  quadrature  leads  to  a  proper  energy  only  for  the 
translational  rigid  body  motion,  and  in  fact  gives  zero  energy 
for  the  uniform  curvature  modes.  For  the  remaining  rigid-body 
and  constant  strain  states,  the  reduced  mass  formulation  underes¬ 
timates  the  kinetic  energy  and  may  predict  frequencies  which  are 
less  accurate  than  the  projection  method. 

Based  upon  the  observations  summarized  in  this  section,  the 
recommended  mass  formulation  for  the  bilinear  Mindlin  plate 
element  therefore  involves  a  single-point  quadrature  for  the 
inplane  motions,  and  the  projection  method  for  the  transverse 
displacements  and  rotations: 

t. 

¥,1]  0  0  0  r2h  ■ 

0  R1H(1)  0  _R2”  0 

Mopt  =  0  0  R1H  0  0  (155) 

0  -r2h  0  r3h  0 

R2H  0  0  0  R3H 

The  additional  effort  required  to  form  the  projected  mass  is 
minimal,  since  the  necessaiy  submatrices  occur  also  in  i-he 
stiffness  calculation.  Furthermore,  the  mass  computation  is 
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usually  performed  only  once  per  analysis,  even  in  nonlinear 
problems,  and  represents  a  negligible  fraction  of  the  complete 
solution  in  all  but  the  smallest  problems. 

The  semi-definite  property  of  the  single-point-integrated 
mass  and  the  projected  mass  matrix  appears  to  present  a  potential 
source  of  difficulty  in  some  methods  of  eigenvalue  extraction, 
such  as  subspace  iteration.  However,  the  subspace  projection  of 
the  mass  will  remain  positive  definite  unless  one  or  more  trial 
vectors  correspond  precisely  to  a  global  deformation  mode  which 
is  free  of  kinetic  energy.  This  situation  is  unlikely  provided 
the  number  of  trial  vectors  is  small  compared  with  the  order  of 
the  system,  which  in  normally  the  case. 
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CHAPTER  6 
MATERIAL  MODELING 


The  scope  of  the  present  investigation  is  limited  to  elastic 
behavior  only.  However,  the  use  of  advanced  composite  materials 
in  turbomachinery  components  is  increasing,  and  effective  methods 
for  analyzing  these  materials  are  needed.  In  this  Chapter  we  in¬ 
troduce  a  technique  for  modeling  multilayered  components  without 
increasing  the  size  of  the  overall  finite  element  model.  The 
approach  leads  to  simple  elements  which  yield  reasonably  accurate 
stress  data,  and  is  applicable  to  most  shear-flexible  structural 
elements. 

6.1  BACKGROUND 

Problems  of  multilayered  plates  and  shells  are  important  in 
.  .  51 

the  design  of  composite  structures,  impact-resistant  vehicle 
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components,  and  vibration-control  treatments;  such  problems 
are  also  of  great  interest  in  the  design  of  the  next  generation 
of  propulsion  system  components.  A  wide  variety  of  theoretical 
and  numerical  treatments  of  such  problems  have  been  developed, 
but  most  of  these  possess  characteristics  which  limit  their  wide¬ 
spread  use  in  production  analysis  software.  As  a  result,  many  of 
the  more  powerful  methods  available  for  analyzing  multilayered 
structures  are  inaccessible  to  analysts  and  designers,  who  must 
resort  to  standard  elements  and  methods  which  are  more  complex 
and  expensive. 

Two  types  of  approaches  predominate  in  the  work  performed  to 

date  in  multilayered  plate  and  shell  analysis.  The  first  of 

these  involves  the  use  of  independent  rotations  or  related 

unknowns  within  individual  layers  to  capture  the  distribution  of 
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transverse  shear  and  normal  stresses  in  detail.  '  This  class 
of  method  works  quite  well  at  the  expense  of  introducing  new 
nodal  variables  whose  number  depends  upon  the  number  of  layers, 
and  which  are  beyond  the  data-handling  scope  of  many  production 
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analysis  codes.  The  second  common  approach  is  through  hybrid 
finite  elements  with  assumed  layer  stress  fields,56-58  with  the 

corresponding  layer  rotations  appearing  as  degrees  of  freedom. 

.  59 

Recently,  Spilker  reported  a  multilayered  hybrid  element  which 
uses  only  six  degrees  of  freedom  per  node,  but  which  is  limited 
to  thin  laminates. 

The  methods  presented  here  deal  effectively  with  multilay¬ 
ered  components,  require  a  minimal  amount  of  added  computation, 
and  may  be  used  in  conjunction  with  most  common  plate  and  shell 
finite  elements.  The  approach  is  based  upon  the  definition  of 
shear  flexibility  corrections  to  be  applied  to  the  basic  plate  or 
shell  element,  and  recovery  of  transverse  shear  stresses  via  the 
equations  of  equilibrium. 

First  we  describe  the  basic  aspects  of  the  shear  flexibility 
correction  as  it  applies  to  layered  isotropic  materials.  We  then 
discuss  the  modifications  which  are  needed  for  some  orthotropic 
laminates,  for  which  a  clear  interpretation  of  the  method  depends 
upon  uncoupling  the  transverse  shear  force  resultants.  Finally, 
procedures  for  point  stress  recovery  are  summarized. 


6.2  LAMINATE  STIFFNESS  CHARACTERISTICS 


The  model  used  herein  is  based  upon  Mindlin's  theory  of 
3  6 

plates.  Section  5.2  outlines  the  kinematic  assumptions  and 
other  pertinent  aspects  of  this  theory.  We  will  work  in  terms  of 


the 
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strains 

and  stresses 
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The  relationship  between  these  generalized  deformation  and  force 
quantities,  as  used  in  Equation  (127) ,  is  often  expressed  in  the 
form: 
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elastic  stiffness  resultants 
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laminated  plate  theory  ;  that 

Aij' 

is : 

Bij, 

and 

Dij  are 

(158) 


defined  as 


(Aij , ) 


t/2  _  2 

Q..(l,z,z  )  dz 
-t/2  30 


(159) 


in  which  are  the  elements  of  the  elasticity  tensor  at  the 
point  in  question,  referred  to  a  common  system  of  coordinates. 


6.3  SHEAR  FLEXIBILITY  CORRECTIONS 


In  most  common  plate  and  shell  elements,  the  assumption  of 
linear  thickness  variations  in  the  tangential  displacements  (see 
equation  110)  results  in  an  extremely  crude  representation  of  the 
transverse  shear  strains.  In  particular,  these  shear  strains  are 
c  istant  through  the  plate  thickness,  and  neither  the  pointwise 
equilibrium  equations  nor  the  traction  boundary  conditions  at  the 
surfaces  are  satisfied  in  general.  For  monolithic,  isotropic 
elements,  a  uniform  reduction  factor  often  is  applied  to  the 
shear  strain  energy  to  obtain  more  realistic  behavior.  Equating 
the  transverse  shear  strain  energy  consistent  with  the  assumed 
displacements  to  that  of  the  parabolic  shear  strain  field  which 
satisfies  the  equilibrium  condition  yields  a  correction  factor  of 
5/6,  which  is  commonly  used  for  isotropic  plates  and  shells. 
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In  the  present  work,  we  use  a  generalization  of  this  idea 

49 

first  proposed  by  Whitney  for  arbitrary  constructions. 

Such  a  correction  is  necessarily  approximate,  but  is  usually 
sufficient  to  bring  the  shear  strain  energy  in  line  with  other 
inodes  of  deformation,  in  a  wav  which  reflects  the  relative  flexi¬ 
bility  of  these  inodes  for  a  given  material  lavup. 

Consider  first  a  layered  construction  for  which  the  shear 
strains  and  resultant  forces  are  related  by: 


Q  " 

k,A,.  0 
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xz 
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1  44 

rxz 
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L  yz  - 
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Based  solely  on  the  elastic  stress-strain  relationship  of  the 
material,  factors  and  k 2  should  both  equal  one.  However,  due 
to  the  excessive  constraint  imposed  by  the  kinematic  assumptions 
of  the  plate  or  shell  theory,  the  strains  7.  produced  by  given 

X  z 

shear  forces  Q^z  are  too  large  over  much  of  the  plate  thickness. 
Accordingly,  the  total  strain  energy  predicted  is  too  large,  and 
the  approximation  appears  too  stiff.  This  error  does  not  respond 
to  mesh  refinement,  since  the  displacement  approximation  through 
the  thickness  remains  linear.  Our  intent  is  to  select  values  for 
k^  and  k2  which  lead  to  stored  energies  of  a  more  reasonable 
magnitude,  and  thus  yield  better  element  behavior. 

Since  the  shear  resultants  are  uncoupled  for  the  case  of  an 
isotropic  material,  the  basic  aspects  of  the  method  can  be  illus¬ 
trated  with  reference  to  a  single  plane.  Below,  we  discuss  the 
determination  of  k^  the  shear  correction  factor  for  the  (x,z) 
plane. 

4  9 

The  shear  corrections  suggested  by  Whitney  depend  upon  the 
assumption  of  cylindrical  bending,  for  which  an  analytical  rela¬ 
tionship  may  be  established  between  the  local  bending  stress  and 
the  transverse  shear  force  resultant: ^ 
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(161) 


a 


(m) 

x,x 


(Bll~Allz)Qxz 


The  superscript  (m)  refers  to  a  particular  layer  within  the 
laminate  cross-section,  and  parameter  D  is  defined  by: 


D 


D11A11 


(162) 


When  combined  with  Equation  (161),  the  equilibrium  equation 


a (m)  +  a(m)  =  o 

x,x  xz,z 


(163) 


can  be  integrated  through  the  plate  thickness  to  obtain  the  shear 
stress  within  a  layer: 
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The  constants  of  integration  a^  are  determined  by  the  condition 

that  a  be  continuous  at  the  layer  interfaces,  and  from  the  free 

surface  boundary  condition  at  either  the  upper  or  lower  surface. 

From  the  condition  that  a  =0  at  z=-t/2,  we  obtain: 

xz 


(1)  -  4SU  Ut+4B11> 


(165) 


in  which  m=l  refers  to  the  bottom  layer  of  the  laminate.  Lettinr 
z^ln^  be  the  lower  surface  of  layer  m,  the  interface  continuity 
conditions  for  m=2,3,...  give: 


(m)  _ 


=  a 
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With  the  above  definitions,  the  strain  energy  density  in  any 
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layer  may  be  written  in  the  form: 

V(m)  =  \  g (m) (2)  Qxz  (167) 

in  which 
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(168) 
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Integrating  Equation  (167)  through  the  laminate  thickness,  and 
equating  the  result  to  the  total  strain  energy  per  unit  area 
obtained  from  Equation  (160) , 
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‘xz 


2klA55 


leads  to  the  shear  correction  factor: 


(169) 


ft/2 

K1  =  ^A44  +  J  g(m)(2)dz]  1  (170) 

The  remaining  factor  k2  may  be  found  in  a  similar  fashion,  using 
the  appropriate  elastic  constants  for  the  (y,z)  plane. 


As  a  representative  example,  consider  a  typical  graphite 
epoxy  material  with  the  properties 


E1 

=  25-xlO6 

E2  =  E3 

=  l.xlO6 

G23 

=  0.2xl06 

G12  =  G13 

=  0.5x10 

V  =  V 

12  13 

=  0.25 

For  this  material,  typical  shear  correction  factors  computed  by 
the  method  outlined  above  are  listed  in  Table  2.  The  classical 
shear  factor  k=5/6  represents  an  upper  bound  under  the  assump¬ 
tions  of  this  Section.  When  small  layers  of  extremely  flexible 
material  are  introduced  in  an  otherwise  uniform  plate,  the  shear 
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Table  2.  Shear  Factors  for  Gr/Ep  Laminates 


Laminate 

kx 

ky 

0.83333 

0.83333 

[0/90] 

0.821 23 

0.821 23 

[ 45/- 45] 

0.68027 

0.68027 

[0/90] 

0.59518 

0.72053 

C-45/45]3 

0.68027 

0.68027 

[-60/0/60] 

0.82579 

0.60800 
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factors  tend  to  drop  quite  rapidly;  this  is  the  case  with  the 
polymeric  materials  used  in  plastic  laminate  interlayers,  and  the 
viscoelastic  materials  typically  employed  in  constrained-layer 
vibration  damping  treatments. 

6.4  UNCOUPLED  CORRECTIONS  FOR  ORTHOTROPIC  LAMINATES 

The  interpretation  of  the  shear  corrections  developed  above 
is  clear  provided  the  transverse  shear  resultants  in  the  (x,z) 
and  (y,z)  planes  are  uncoupled.  When  the  transverse  shear  moduli 
in  these  two  planes  differ,  and  when  layers  with  orientations 
other  than  0°  and  90°  are  present,  the  stress-strain  relation  has 
the  form: 
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which  we  will  represent  by  Q  =  Sy.  For  such  cases,  the  interpre¬ 
tation  of  the  correction  factor  derived  from  cylindrical  bending 
assumptions  is  open  to  question. 


The  existence  of  a  positive  definite  strain  energy  function 
implies  that  the  shear  stiffness  matrix  S  is  real,  symmetric,  and 
positive  definite.  Therefore,  there  exists  a  planar  transforma¬ 
tion  of  coordinates  defining  an  alternate  system  of  reference 
axes  (x',y'),  for  which  the  corresponding  stiffness  S#  is  diag¬ 
onal  and  the  shear  resultants  are  uncoupled.  Since  the  z  axis 
remains  unchanged,  the  transformation  by  an  angle  /)  has  the  form: 


vith 


S'  =  qsqt 


cos0  -sin$ 
sin/3  cos/J 


(172) 

(173) 


The  required  angle  of  rotation  is  easily  determined  in  terms  of 
the  original  shear  stiffness  coefficients: 
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tan  (2/5) 
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(174) 


6.5  SHEAR  STRESS  RECOVERY 


For  pointwise  stress  recovery  consistent  with  the  transverse 
shear  flexibility  correction  outlined  here,  Equation  (164)  may  be 
used  directly  at  any  station  z  within  the  laminate  thickness, 
with  the  shear  force  resultants  obtained  from  Equation  (160). 

The  integration  constants  a^  may  be  stored  and  recalled  for  use 
in  the  stress  recovery,  at  a  cost  of  only  one  floating  point  word 
per  layer.  With  the  assumption  of  linear  displacement  variation 
through  the  plate  thickness,  the  predicted  transverse  shear 
stress  field  is  quadratic  within  each  layer. 

For  the  bilinear  plate  element  used  in  the  present  work,  we 
choose  to  evaluate  the  transverse  shear  stresses  at  the  element 
center,  which  corresponds  to  the  optimal  sampling  point.  With 
higher-order  displacement  elements,  it  may  be  possible  to  obtain 
accurate  transverse  shear  stresses  at  a  regular  grid  of  points 
within  an  element,  such  as  the  2x2  Gauss  points.  These  data  may 
be  used  in  turn  for  the  evaluation  of  transverse  normal  stresses 
on  a  smaller  grid,  by  integration  of  the  remaining  equilibrium 
equation. 
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CHAPTER  7 

NUMERICAL  EXAMPLES 


This  Chapter  presents  a  number  of  solved  problems  which 
demonstrate  the  analytical  methods  discussed  earlier.  Many  of 
these  are  small,  relatively  simple  problems  for  which  closed  form 
solutions  or  previous  numerical  results  exist.  Comparisons  with 
known  results  are  made  both  to  verify  specific  capabilities  and 
to  determine  the  accuracy  characteristics  of  the  present  analysis 
techniques . 

The  numerical  solutions  are  presented  in  four  sections, 
which  pertain  to  four  primary  areas  of  investigation  in  the  work 
performed.  Section  7.1  deals  with  basic  dynamic  problems,  and 
illustrates  the  performance  of  the  stabilization  procedure  and 
the  optimal  mass  formulation  for  the  bilinear  plate  element  (see 
Chapter  5).  In  Section  7.2,  we  present  several  problems  involv¬ 
ing  composite  materials  or  layered  wall  construction;  all  of 
these  are  solved  using  a  single  layer  of  plate  elements,  based  on 
the  shear  correction  technique  described  in  Chapter  6.  Section 
7.3  contains  a  number  of  sensitivity  analyses  (Chapter  3),  and 
Section  7.4  presents  analyses  using  the  probabilistic  techniques 
of  Chapter  4. 

7 . 1  DYNAMICS  EXAMPLES 


The  problems  of  this  Section  demonstrate  the  bilinear  plate 
element  formulation  of  Chapter  5.  In  particular,  we  contrast  the 
recommended  mass  formulation  with  other  commonly-used  alternative 
forms.  The  first  example  is  an  axial  vibration  problem  which  is 
simple  in  concept,  but  which  illustrates  very  well  the  ability  of 
the  present  mass  formulation  to  move  spurious  dynamic  modes  to 
the  top  of  the  frequency  spectrum.  The  remaining  two  problems 
are  cases  for  which  troublesome  results  have  been  reported  in  the 
past;  the  present  method  gives  a  reliable  solution  with  improved 
accuracy. 
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7.1.1  Comparison  of  Mass  Formulations  for  Axial  Vibration 

This  example  demonstrates  the  occurrence  of  artificial  modes 
controlled  by  the  anti-hourglass  stiffness  parameters.  Consider 
the  planar  vibrations  of  a  long,  thin  strip  as  shown  in  Figure 
12.  One  quadrant,  with  dimensions  (1,10,0.05),  is  modeled  by  a 
single  element;  symmetry  is  imposed  along  the  axes  x=0  and  y=0. 

.  7 

The  mechanical  properties  are  E=10  ,  i/=0.25,  and  p=0. 000259. 

Solutions  have  been  calculated  with  all  four  of  the  mass 
formulations  discussed  previously,  and  results  are  listed  in 
Table  3.  Predicted  frequencies  corresponding  to  the  inplane 
staircase  mode  are  shown  in  parentheses.  Accurate  estimates  of 
the  exact  natural  frequencies  are  not  expected,  due  to  the  coarse 
mesh  employed;  our  intent  in  this  example  is  to  examine  the 
occurrence  of  spurious  low-frequency  modes  for  each  of  the  mass 
formulations. 

The  solution  with  lumped  mass  exhibits  the  lowest  spurious 
frequencies,  as  well  as  extreme  lower  bounds  on  the  first  two 
physical  modes,  since  half  of  the  element  mass  is  concentrated  at 
the  free  end.  The  consistent  and  projected  masses  give  similar 
frequency  estimates,  with  a  first  spurious  mode  quite  close  to 
the  real  fundamental  frequency;  both  spurious  modes  in  these  two 
solutions  are  inplane  staircase  modes,  and  the  slightly  lower 
frequencies  predicted  with  full  consistent  masses  are  due  to  the 
nonzero  hourglassing  kinetic  energy. 

The  single-point  integrated  mass  yields  a  reliable  solution, 
though  the  natural  frequencies  are  somewhat  higher  than  with  the 
consistent  mass  and  the  projection  method.  The  mass  for  this 
case  has  been  augmented  by  a  small  fraction  (0.1%)  of  the  lumped 
mass  to  achieve  positive  definiteness,  since  all  the  modes  a:  ■». 
solved.  The  two  highest  frequencies  are  controlled  by  the  lumped 
mass  contribution. 
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Figure  12.  Slender  Strip  Geometry  and  Properties. 
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Table  3.  Comparison  of  Results  for  Planar  Vibration  of  Thin  Strip 


Mass  (1*2  w3  u>4 

Lumped  (16,757.)  27, 779. a  (177,416.)  287 ,083. b 

Consistent  34,023. a  (35,547.)  351, 611. b  (376,364.) 

Projected  34,023. 3  (41,046.)  351, 611. b  (434,582.) 

Reduced  39, 277. 3  405, 904. b  (1,059,812.)  (11,220,998.) 


3  Stretching 
b  Stretching 


mode,  long  direction, 
mode,  short  direction. 
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Table  3  shows  an  additional  solution  using  the  recommended 
mass  technique  (lxl  quadrature)  and  four  elements.  The  first 
stretching  mode  is  quick  to  converge,  suggesting  that  the  poor 
one-element  solution  is  not  indicative  of  an  unforeseen  pathology 
in  the  element. 

The  energy  content  of  the  mode  shapes  for  this  example  is 
unambiguous,  and  the  non-physical  vibration  modes  could  have  been 
rejected  automatically  on  this  basis.  However,  the  occurrence  of 
such  spurious  solutions  in  larger  problems  may  limit  the  number 
of  legitimate  modes  obtained,  and  exacts  added  storage  demands 
which  may  be  unacceptably  large. 

7.1.2  Vibration  of  a  Corner-Supported  Plate 

The  vibrations  of  a  corner-supported  plate  (Figure  13)  have 
been  considered  by  Belytschko  and  Tsay45  using  the  stabilized 
Mindlin  plate  element.  Reference  45  contains  results  for  the 
first  three  frequencies,  as  functions  of  the  w  and  6  hourglass 
stiffnesses;  artificial  or  inaccurate  frequencies  were  obtained 
only  when  one  or  both  of  these  stiffnesses  were  suppressed. 

In  our  analysis,  we  assume  double  symmetry  and  consider  out- 
of-plane  motions  only.  The  length  of  each  edge  of  the  entire 
plate  is  24;  the  material  properties  are  E=360,000,  y=0.38,  and 
p=0 . 001 .  A  uniform  mesh  of  36  elements  is  used,  so  that  natural 
frequencies  should  be  directly  comparable  with  those  of  Reference 
45. 


Table  5  lists  normalized  frequency  values  for  the  first  fivr 
symmetric  vibration  modes  of  the  plate.  All  values  are  in  reas¬ 
onable  agreement  with  the  exact  results,  and  with  the  numerical 
values  predicted  by  Belytschko  and  Tsay;  the  minor  differences 
which  do  exist  in  the  numerical  solutions  can  be  attributed  to 
the  differing  parameters  used  in  the  stiffness  stabilization. 
Recall  that  the  recommended  technique  for  bending  motions  uses 
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Table  4.  Vibration  Modes  of  Thin  Strip  (Four-Element  Solution) 


Mode 

Predicted 

Frequency 

Exact 

Frequency 

Description  of 

Mode 

1 

31  ,259. 

30,865. 

y 

stretching, 

first 

mode 

2 

104,840. 

92,596. 

y 

stretching, 

second 

1  mode 

3 

232,755. 

154,326. 

y 

stretch ing , 

third 

mode 

4 

405,266. 

308653- 

X 

s  tretch ing, 

first 

mode 
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E  ■  430,000 
v  *  0.38 

p  =  0.001 

t  =  0.375 


Table  5.  Natural  Frequencies  for  Corner-Supported  Plate. 

a)  -  wa2(D/pt)  ^ ^ 


Mode 

Consi stent 

Mass 

Projection 

Method 

Reduced 

Quadrature 

Bely tschko 

&  Tsay  [45] 

Analytic 

1 

7.118 

7.118 

7.124 

7.099-7.185 

7 . 1  20 

2 

18.79 

1  8.79 

19.08 

19.18  -19.19 

19.60 

3 

44.01 

44.01 

44.79 

42.70  -43.98 

44.40 

4 

95.18 

95.33 

98.11 

- 

- 

5 

124.13 

124.14 

1 32.44 

“ 

- 
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the  projection  method.  For  this  case,  the  mass  matrix  obtained 
by  single  point  quadrature  leads  to  slightly  higher  frequencies, 
due  to  neglect  of  the  kinetic  energy  associated  with  constant- 
curvature  states. 

7.1.3  Vibration  of  a  Free-Free  Square  Plate 

This  example,  also  taken  from  Reference  45,  exhibits  a  free 
vibration  mode  which  is  sensitive  to  the  w-hourglass  stiffness. 
The  geometry  and  properties  are  identical  to  those  in  the  pre¬ 
vious  example,  but  the  entire  plate  is  modeled.  A  36-element 

mesh  is  used,  to  facilitate  comparisons  with  the  solutions 

4  5 

reported  by  Belytschko  and  Tsay. 

Table  6  summarizes  the  normalized  frequencies  obtained  for 
the  first  six  bending  modes  of  the  plate  (the  first  three  modes, 
which  correspond  to  rigid-body  motions,  are  not  listed) .  The 
projection  method  is  clearly  superior,  particularly  for  the  third 
frequency  which,  according  to  Reference  45,  is  quite  sensitive  to 
the  hourglass  stiffness  parameter.  This  frequency  corresponds  to 
the  (3,1)  mode  of  the  plate;  the  accuracy  is  particularly  good  in 
view  of  the  fact  that  three  half-waves  are  represented  by  only 
six  bilinear  elements. 

7.2  COMPOSITES  AND  LAYERED  STRUCTURES 

The  examples  of  this  Section  involve  both  advanced  composite 
materials  and  layered  (sandwich)  components.  They  illustrate  the 
use  of  the  shear  corrections  of  Chapter  6  for  static  and  dynamic 
problems.  The  overall  stiffness  effect  is  very  realistic,  and 
good  results  can  be  expected  for  resultant  (in  the  sense  of  plate 
theory)  quantities.  Two  examples  present  point  stress  results; 
these  are  of  reasonable  quality,  but  seem  somewhat  more  vul¬ 
nerable  to  the  errors  entailed  in  the  assumption  of  cylindrical 
bending  than  the  stiffness  properties. 


Table  6.  Natural  Frequencies  for  Free-Free  Plate. 

w  ■  uia2  ( D/pt )  ^ 


Mode 

Projection 

Method 

Reduced 

Quadrature 

Bely tschko 

&  Tsay  [45] 

Analytic 

1 

(22) 

13.07 

13.42 

13.14 

13.47 

2 

(13) 

19.14 

20.46 

18.12 

19  .60 

3 

(31  ) 

25.81 

27.46 

19.05 

24.27 

4 

(32) 

34.11 

36.85 

- 

35.02 

5 

(23) 

34.11 

36.85 

- 

35.02 

6 

(41  ) 

62.87 

70.59 

- 

61.53 

7.2.1  Unsyraetric  Laminated  Plate 

The  semi-infinite  thick  plate  shown  in  Figure  14  is  sub¬ 
jected  to  a  sinusoidal  pressure  load  q(x)  =  qQsin(7tx/a)  ,  and  is 
simply  supported  on  its  lateral  edges.  The  0*  direction  is  the 
fiber  direction  in  the  top  layer,  and  corresponds  to  the  infinite 
(y)  direction.  The  material  properties  are  EL/ET=25,  GTT/ET=0.5, 
Gtt/Et=0.2,  and  v TT=0 *25«  The  plate  has  a  width  2a=24  and 
thickness  t=6.  An  exact  elasticity  solution  of  this  problem  has 
been  presented  by  Pagano;  finite  element  results  based  upon  the 

use  of  independent  layer  rotations  are  reported  by  Palazotto  and 
54 

Witt. 


The  plate  is  modeled  using  ten  elements  over  half  the  width, 
with  symmetry  conditions  applied  at  the  centerline.  Transverse 
shear  stresses  in  the  element  nearest  the  support  are  shown  in 
Figure  15.  The  results  are  in  reasonable  agreement  with  the 
exact  solution,  with  the  peak  shear  stress  being  overestimated  by 
about  eight  percent.  The  finite  element  solution  of  Reference 
54,  using  30  elements  with  independent  rotations  in  each  layer, 
appears  to  overestimate  the  maximum  shear  stress  by  three  to  four 
percent,  based  on  graphical  results  presented  therein. 

7.2.2  Three-Layered  Plate  under  Pressure 

The  square  plate  in  Figure  16  is  a  [0/90/0]  graphite/epoxy 

laminate,  with  EL=25xio6  and  the  remaining  properties  defined  as 

in  the  previous  example.  The  0°  direction  is  aligned  with  the  x 

axis.  For  the  sinusoidal  pressure  q(x,y)  =  q.sin  (flx/a)  sin  (7ry/b)  , 

^  .62 

and  simply  supported  edges,  an  analytical  solution  is  possible. 

We  consider  the  case  a=b=10,  for  plate  thicknesses  between  0.1 
and  2.5.  In  the  finite  element  model,  an  llxii  mesh  of  bilinear 
elements  is  used  to  represent  the  entire  plate.  The  symmetry  of 
the  problem  is  not  exploited,  in  the  interest  of  obtaining  stress 
values  at  suitable  locations  for  comparison  with  other  solutions. 
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NORMALIZED  SHEAR  STRESS,  «xz/Qo 


Figure  15.  Transverse  Shear  Stresses  in  Unsymmetric  Plate 
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Figure  16.  Square  [0/90/0]  Plate  under  Pressure  Load. 
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Table  7  compares  the  normalized  bending  stresses  (defined  by 

—  22.. 

o=o t  /qna  )  obtained  using  the  present  method  with  the  elasticity 
solution  and  the  finite  element  results  of  Engblom  and  Ochoa, 
who  use  a  40-DOF  plate  element  with  higher-order  displacement 
variations  through  the  thickness.  Results  listed  in  the  Table 
are  limited  to  those  values  reported  in  Reference  62  which  can  be 
evaluated  directly  at  element  centers. 

The  bending  stresses  obtained  from  the  present  solution  are 
in  reasonable  agreement  with  the  remaining  solutions.  The  trends 
predicted  are  quite  similar  to  those  of  the  higher-order  element, 
but  the  accuracy  obtained  is  generally  lower.  In  this  example, 
the  deformation  pattern  is  quite  different  from  the  cylindrical 
bending  assumptions  used  to  derive  the  shear  corrections.  As  a 
result  the  deflections  and  overall  load  paths  are  reaonable,  but 
pointwise  stress  accuracy  is  limited. 

7.2.3  Circular  Sandwich  Plate 

The  finite  element  mesh  in  Figure  17  represents  one  quadrant 
of  a  circular  sandwich  panel  with  the  following  properties: 

Faces:  E  =  lxlO7  v  =  0.30  tf  =  0.025 

Core:  G  =  260.000.  t  =  0.450 

-  c 

The  radius  of  the  panel  is  a=20,  and  the  outer  edge  is  completely 

fixed.  A  uniform  static  pressure  q  =  10  is  applied.  This  case 

64 

has  been  analyzed  by  Sharifi  ,  using  special-purpose  elements 
with  independent  shear  rotations  in  the  sandwich  core  layer. 

Figure  18  shows  the  radial  distribution  of  moment  resultants 
obtained  from  the  present  analysis.  Though  only  graphical  re¬ 
sults  are  available  in  Reference  64,  the  two  solutions  appear  to 
agree  quite  well.  The  shear  force  per  unit  length  obtained  from 
the  finite  element  solution  is  shown  in  Figure  19,  together  with 
the  exact  solution  Q(r)  =  qr/2. 
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Table  7.  Normalized  Stresses  for  Square  [0/90/0]  Plate 


Stress 

Elastic ity 

Higher-Order 

a/h 

Component  Present  Solution  [62] 

Element  [63] 

4 

_  1 

°x 

0.357 

0.755 

0.391 

°y 

0.521 

0.556 

0.572 

10 

°x 

0.476 

0.590 

0.500 

°y 

0.273 

0.285 

0.279 

20 

°x 

0.506 

0.552 

0.531 

>> 

1  o 

0.200 

0.189 

0.210 

50 

°x 

0.513 

0.541 

0.541 

°y 

0.176 

0.185 

0.164 

100 

°x 

0.514 

0.539 

0.542 

a 

°y 

0.173 

0.181 

0.167 

'  5x  at 

center  of 

plate,  z»-| 

2  oy  at 

center  of 

plate  ,  z-|-  ( top 

of  90°  layer) 
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lqure  17.  Circular  Sandwich  Plate. 


MOMENT  RESULTANTS 


Figure  18.  Moment  Resultants  in  Circular  Sandwich  Plate. 
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Figure  19.  Shear  Forces  in  Circular  Sandwich  Plate. 
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7.2.4  Rectangular  Sandwich  Plate 

A  square  sandwich  panel  (Figure  20)  is  loaded  by  a  uniform 

pressure  qQ.  The  three-layer  plate  is  50  inches  on  each  side, 

with  identical  aluminum  face  sheets  (E=10.5xio6,  i/=0.3,  tf=0.015) 

and  a  honeycomb  core  (G=50,000,  t  =1.0).  All  edges  of  the  panel 

c 

are  completely  fixed. 

The  present  solution  uses  a  5x5  mesh  of  bilinear  elements  in 
one  quadrant  of  the  plate,  and  yields  a  transverse  deflection  at 

the  center  w  =  0.09285.  This  value  compares  well  (2.2  percent) 

^  65 

with  the  analytic  solution  presented  by  Kan  and  Huang,  which 

gives  wc  =  0.09497.  The  finite  element  solutions  of  References 

66  and  67  achieve  comparable  accuracy,  but  with  more  than  twice 

as  many  equations  and  considerably  more  complicated  elements. 

7.2.5  Vibrations  of  a  Layered  Panel 

The  natural  frequencies  of  a  rectangular  [0/90/0]  laminate 

obtained  using  the  present  analysis  have  been  compared  with  the 

.  .  34 

analytical  solution  by  Ashton  and  Whitney.  The  plate  (Figure 

21)  has  dimensions  30xio,  and  mechanical  properties  identical  to 
those  used  in  the  first  two  examples v  a  density  of  p  =  0.0001  is 
assumed.  Each  layer  is  0.01  thick.  In  the  finite  element  solu¬ 
tion,  a  15x5  element  mesh  is  used  to  represent  the  entire  plate. 
All  four  edges  are  simply  supported. 


Table  8  compares  the  computed  natural  frequencies  with  the 
analytical  solution, 


D 


5  -  “V'Vd22  ‘  1 5^) 


llr-b^  .  2  _(D.12.+  2.D66)  r  ^bi2  A  ^4,1/2 


D 


22 


-[ran;)*  +  n  ] 


(175) 


in  which  (m,n)  are  mode  numbers  along  the  (x,y)  axes.  The  seven 
lowest  frequencies  in  the  Table  represent  all  modes  below  the 
first  n=3  mode.  The  n=l  modes  exhibit  good  accuracy;  for  n=2 , 
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Figure  20.  Clamped  Sandwich  Panel  under  Uniform  Pressure. 
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Table  8.  Natural  Frequencies  of  [0/90/0]  Plate 


* 

Mode 

m 

n 

^exact 

U) 

comp. 

Error(  51) 

1 

1 

1 

1.415 

1.473 

4.1 

2 

2 

1 

2.626 

2.733 

4.1 

3 

1 

2 

4.420 

4.864 

10.0 

4 

3 

1 

4.622 

5.056 

9.4 

5 

2 

2 

5.659 

6.358 

12.4 

6 

4 

1 

7.406 

8.025 

8.4 

7 

3 

2 

7.691 

8.528 

10.9 
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where  the  five  elements  across  the  width  can  be  expected  to  give 
only  marginal  accuracy,  the  computed  frequencies  are  still  within 
10-12  percent  of  the  exact  values. 

7.3  SENSITIVITY  ANALYSIS  EXAMPLES 

The  examples  of  this  Section  are  sensitivity  calculations, 
in  which  material  modulus  and  density,  plate  thickness,  and 
arbitrary  geometric  variables  appear  as  independent  parameters. 
The  first  five  problems  have  analytical  solutions,  so  that  errors 
in  the  finite  element  solution  can  be  assessed  conclusively.  In 
these  problems,  we  include  both  static  and  natural  frequency 
results,  and  sensitivities  with  respect  to  intrinsic  properties, 
shape  (dimensions) ,  and  orientation.  The  final  example  deals 
with  a  twisted  plate  for  which  numerical  results  are  available, 
and  compares  the  sensitivity  calculations  for  total  angle  of 
twist  to  approximations  obtained  using  finite  differences. 

7.3.1  Static  Analysis  of  a  Tension  Strip 

The  long,  thin  strip  in  Figure  12  (see  Section  7.1.1)  is 
subjected  to  a  uniform  load  applied  at  the  end.  We  choose  as 
sensitivity  variables  the  modulus  E,  thickness  t,  width  b,  and 
length  L.  The  first  two  of  these  are  intrinsic  variables,  and 
the  remaining  two  are  geometry  parameters  which  affect  the  nodal 
positions.  The  exact  inplane  displacements  are  linear  functions 
of  position,  as  are  the  sensitivities,  and  therefore  a  single 
bilinear  element  should  reproduce  both  results  exactly.  Data 
obtained  for  the  displacement  and  stress  resultant  sensitivities 
from  a  single-element  model  are  indeed  exact,  as  shown  in  Table 
9. 

7.3.2  Statics  of  a  Cantilever  Beam 

Figure  22  shows  a  cantilever  beam  subjected  to  a  transverse 
force  at  the  tip.  Again,  an  analytical  solution  is  possible  both 
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Table  9.  Sensitivity  Data  for  Simple  Tension  Problem 


Quantity 

Exact 

Exact  (x-L) 

Computed 

u 

Px/Ebt 

0.001 

0.001 

3u/3E 

-Px/E2bt 

i 

o 

X 

O 

o 

1 

• 

o 

X 

o 

1 

o 

3u/3t 

-Px/Ebt2 

-0.01 

-0.01 

3u/3b 

-Px/Eb2t 

-0.001 

-0.001 

3u/3L 

P/Ebt 

0.0001 

0.0001 

N 

P/b 

100.0 

100.0 

3N/3E 

0 

o 

• 

o 

o 

o 

3N/3t 

0 

o 

• 

o 

o 

o 

3N/3b 

-P/b2 

-100.0 

-100.0 

3N/3L 

0 

0.0 

0.0 

1 07 ;  v-0.3; 

t-0.1;  L - 1 0  ; 

b-1;  P-100 

Figure  22.  Cantilever  Beam  with  Tip  Load 
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for  the  displacement  and  rotation,  and  for  the  sensitivities  with 
respect  to  modulus  E,  thickness  t,  and  width  b.  Five  bilinear 
plate  elements  are  used  to  model  the  beam;  this  number  is  suffi¬ 
cient  for  good  accuracy  but,  since  the  elements  have  only  linear 
displacement  and  rotation  fields,  does  not  reproduce  the  exact 
solution.  Table  10  summarizes  computed  results  for  the  displace¬ 
ments  and  rotations.  Note  that  the  displacement  sensitivities 
are  no  less  accurate  than  the  displacements  themselves  (all  are 
approximately  1  percent  in  error) ,  and  that  the  rotational 
results  are  exact.  Table  11  shows  moment  and  shear  results,  and 
the  force  sensitivities  which  are  nonzero.  In  all  cases,  the 
moment  and  shear  sensitivities  are  exact,  despite  small  errors  in 
the  displacement  solution.  It  is  probably  reasonable  to  expect 
exact  results  for  the  force  sensitivities  in  a  statically 
determinate  problem. 

7.3.3  Orientation  Sensitivity  of  a  Beam 

Consider  the  bar  shown  in  Figure  23,  which  is  inclined  with 
respect  to  the  global  X  axis  and  subjected  to  a  vertical  force, 
producing  both  stretching  and  bending  response.  This  problem  is 
intended  to  test  the  sensitivity  calculation  for  element  local 
axis  orientation;  we  select  the  orientation  6  as  the  geometric 
control  variable,  which  leads  to  B'=0,  | j| '=0,  and  A'*0  for  all 
elements.  Note  that  the  nodal  coordinate  sensitivities  are 
simply  X ' =-Y ,  Y'=X. 

For  S-0,  the  nonzero  results  and  sensitivities  are  given  in 
Table  12.  Computed  results  are  obtained  from  a  model  with  five 
bilinear  Mindlin  plate  elements,  which  exhibits  a  moderately 
small  displacement  error  (1-3  percent).  Again,  the  displacement 
sensitivities  exhibit  errors  which  are  similar  in  magnitude  to 
the  displacement  error  in  the  original  solution;  both  the 
computed  moments  and  axial  force  sensitivities  are  exact. 
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Table  10. 


Displacement  Sensitivity  Data  for  Cantilever  Beam 


Quantity 

Exact 

Exact  (x-L) 

Computed 

w 

2P(3Lx2-x3) /Ebt3 

0.4 

0.39602 

3w/3E 

2  7,  2  2. 

-2P(3Lx-xJ)/Ebt"5 

-4.0*10~8 

-3.96k10-8 

3w/3t 

-6P(3Lx2-x3)/Ebt4 

-1  2.0 

-11.88 

3w/3b 

-2P(3Lx2-x3) /Eb2t3 

-0.4 

-0.39602 

e 

6P (2Lx-x2) /Ebt3 

0.06 

0.06 

38/3E 

-6P(2Lx-x2)/E2bt3 

-6.0*10~9 

-6. 0*1 0~9 

39  /  3t 

-1 8P(2Lx-x2)/Ebt4 

-1  .8 

-1  .8 

3  0  /3b 

-6P(2Lx-x2)/Eb2t3 

-0.06 

-0.06 

1-1*10;  v 

-0;  t-0.1;  b-1;  L - 1  0  ; 

P-1 
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Table  11.  Force  Sensitivity  Data  for  Cantilever  Bean/ 


Element  Centers 


Quantity 

x-1 

x-  3 

x-5 

x- 7 

x-9 

M  Exact 

-9.0 

-7.0 

-5.0 

-3.0 

-1  .0 

Comp. 

-9.0 

1 

— -3 

• 

O 

-5.0 

-3.0 

-1  .0 

3M/3b  Exact 

9.0 

7.0 

5.0 

3.0 

1  .0 

Comp . 

9.0 

7.0 

5.0 

3.0 

1  .0 

Q  Exact 

1  .0 

1 .0 

1  .0 

1.0 

1  .0 

Comp . 

1  .0 

1.0 

1.0 

1.0 

1  .0 

3Q/3b  Exact 

-1  .0 

-1  .0 

-1  .0 

-1  .0 

-1  .0 

Comp. 

-1  .0 

-1  .0 

-1  .0 

-1 .0 

-1.0 

E-1*107;  v-0;  t-0.1; 

b-1; 

L - 1 0 ;  P-1 

3Q 

'  3S  * 

3Q 

3 1 

Only  nonzero  values 

shown 

3M  3M 

;  TE  *  3t 

-  0 

identically 

Table  12.  Results  for  Angular  Orientation  Problem  (9-0)* 


Quantity 

Nodal  Positions 

x-L/5 

X-2L/5 

X-3L/5 

X-4L/5 

x-L 

v  Exact 

0.07000 

0.26000 

0.54000 

0.88000 

1 . 2500 

Comp . 

0.06756 

0.25512 

0.53268 

0.87024 

1 . 2378 

au/38  Exact 

-0.06998 

-0.25995 

-0.53993 

-0.87990 

-1 . 2499 

Comp. 

-0.06754 

-0.25507 

-0.53260 

-0.87014 

-1 .2377 

Element  Centers 

Quantity 

x-L/10 

X-3L/10 

x-L/2 

X-7L/10 

X-9L/10 

M  Exact 

-4.500 

-3.500 

-2.500 

-1 .500 

-0.500 

M  Comp. 

-4.500 

-3.500 

-2.500 

-1 .500 

-0.500 

3N/30  Exact 

1 .000 

1 .000 

1  .000 

1  .000 

1 .000 

3N/39  Comp. 

1 .000 

1 .000 

1 .000 

1  .000 

1 .000 

*E«400 , 000 ;  v-0;  t-0.1;  b-1;  L-5;  F-1. 
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Table  13  summarizes  the  results  for  a  similar  calculation  with 
0=26 . 565*  (tan0=l/2) .  Computed  and  exact  displacement  values 
again  compare  well.  Results  for  the  axial  force  and  moment 
resultants  are  exact  as  above  (only  one  set  of  values  is  shown  in 
the  Table) . 

7.3.4  Frequency  Sensitivity  of  a  Flat  Strip 

The  planar  strip  of  Figure  12  is  considered  again,  to 
determine  parameter  sensitivities  of  the  fundamental  frequency. 
Using  a  single  bilinear  element  (which  provides  a  poor  estimate 
of  the  lowest  frequency) ,  we  can  make  some  interesting  observa¬ 
tions  on  the  sensitivity  solution.  Since  the  first  natural 

2  2 

frequency  is  u-jEn  /4/>L  ,  we  note  that  dw/3E=w/2E,  du/3t=0,  and 

3u/3  Lr=-u  /L* 

7 

The  lowest  natural  frequency  for  a  strip  with  E=10  ,  i/=0.3, 
L=10,  and  p=0 . 000259,  as  computed  using  a  single  element  with 
consistent  mass,  is  uc=39, 669.4  (see  Section  7.1.1),  and  compares 
poorly  with  the  exact  value  of  w=30,865.3.  The  sensitivities, 
compared  with  exact  results,  are  similarly  poor.  However, 
sensitivity  values  computed  on  the  basis  of  the  finite  element 
model  frequency  (e.g.,  du/dE=u /2E)  are  nearly  exact,  as  shown  in 
Table  14.  That  is,  the  sensitivity  values  are  related  to  the 
model  frequency  in  the  correct  manner,  and  the  errors  in  the 
computed  sensitivities  are  dominated  by  the  discret ization  error 
in  the  original  solution.  This  is  true  because  parameters  such 
as  the  modulus  and  density  (and,  in  this  problem,  the  length) 
enter  the  finite  element  solution  in  precisely  the  same  way  as 
for  the  analytical  problem. 
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0  t 

Table  13.  Results  for  Angular  Orientation  Problem  (0-26.565  ) 


Nodal  Positions 

Quantity 

x-L/5 

X-2L/5 

X-3L/5 

X-4L/5 

x-L 

u  Exact 

-0.03912 

-0.14532 

-0.30184 

0.491 89 

-0.69872 

Comp. 

-0.03775 

-0.1 4258 

-0.29772 

0.48641 

-0.69186 

v  Exact 

0.07338 

0.27254 

0.56603 

0.92241 

1  .3102 

Comp . 

0.07553 

0.28552 

0.59553 

0.97293 

1 . 3839 

3u/30  Exact 

-0.05868 

-0.21 798 

-0.45275 

0.73784 

-1 .0480 

Comp . 

-0.05662 

-0.21 387 

-0.44659 

0.72961 

-1 .0378 

3v/36  Exact 

-0.07824 

-0.29064 

-0.60367 

0.98378 

-1 .3974 

Comp . 

-0.07550 

-0.28516 

-0.59545 

0.97281 

-1 . 3837 

Element  Centers 

Quantity 

x-L/10 

X-3L/10 

x-L/2 

X-7L/10 

X-9L/10 

N 

0.4472 

0.4472 

0.4472 

0.4472 

0.4472 

M 

-4.500 

-3.500 

-2.500 

-1 .500 

-0.500 

3N/30 

0.8944 

0.8944 

0.8944 

0.8944 

0.8944 

3M/30 

2.250 

1  .750 

1  .  250 

0.750 

0.250 

f  E- 400 , 000  ; 

o 

■ 

4-> 

o 

1 

> 

;  b  »  1  ;  L  ■  5  . 

5901  7;  F  -1. 

Only  computed  element 

results  are  shown;  all 

values  are 

exact . 
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Table  14.  Frequency  Sensitivities  for  Axial  Vibration  Problem 


2 


Quantity 

Exact 

Exact  Value 

Computed 

3u>/3E 

u>/2E 

0.0019834 

0.0019835 

A 

3u/3t 

0 

0. 

-1 .81 *10'y 

3u/3  L 

-o j/L 

-3,966.94 

-3,965.05 

‘e-UIO7;  v-0.3;  t-0.1;  p-0.000259;  L-10;  b-1. 

2 

Exact  sensitivities  computed  using  F.  E.  model  frequency. 
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7.3.S  Frequency  Sensitivity  of  a  Bean 


For  the  cantilever  beam  (Figure  22) ,  an  analytical  solution 
for  parameter  sensitivities  of  the  natural  frequencies  is  quite 
simple.  Defining 


(176) 


The  bending  frequencies  are  u ^=a ^0 ,  where  are  independent  of 

the  geometry  and  properties  of  the  beam.  In  particular,  the 

first  three  natural  frequencies  have  a  -  3.52,  22.0,  and  61.7, 

,  68 

respectively. 


Table  15  summarizes  the  results  obtained  for  a  particular 
case,  using  three  different  meshes.  The  sensitivity  parameters 
are  modulus  E,  density  p,  thickness  t,  width  b,  and  length  L.  It 
is  instructive  to  study  the  results  from  a  relatively  coarse 
model  (five  bilinear  elements)  first;  this  model  is  labeled  Mesh 
1  in  the  Table.  All  computed  results  for  the  first  mode  are 
quite  good  for  Mesh  1,  with  the  error  in  frequency  sensitivities 
being  similar  in  magnitude  to  the  frequency  error  itself.  For 
the  next  two  modes,  the  sensitivities  for  intrinsic  parameters  E, 
p,  and  t  are  at  least  equal  in  accuracy  to  the  frequencies,  for 
reasons  explained  in  the  last  example.  The  length  sensitivity  in 
Mesh  1  has  been  defined  by  attributing  coordinate  sensitivities 
only  to  the  end  nodes,  however,  and  is  rather  poor:  this  "local" 
geometry  parameter  does  not  enter  the  finite  element  frequency 
equation  in  the  same  manner  as  in  the  analytical  solution. 


Meshes  2  and  3  represent  the  two  obvious  solutions  to  the 
poor  accuracy  of  Mesh  1  for  the  length  parameter  in  higher  modes. 
Mesh  2  is  a  refined  model,  in  which  ten  elements  are  used,  and 
the  coordinate  sensitivities  are  defined  for  the  end  nodes  only, 
as  in  Mesh  1.  All  results  are  much  improved,  as  expected;  but 
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Table  15. 


Frequency  Sensitivities  for  Cantilever  Beam 


t 


Mode 

ID 

3ti)/3E 

3w/3  p 

3u>/3t 

3w/3b 

3u/3L 

1 

Exact 

201 .6 

1 .01 *10~5 

-3.97-105 

2016. 

0. 

-40.32 

Mesh-1 

202.3 

1 .01 *io'5 

-3.98*105 

2023. 

3.9*10' 

10 

-40.67 

Mesh -2 

201.6 

1 .01 *io'5 

-3.97*105 

2016. 

4 .9  *1  o' 

10 

-40.37 

Mesh-3 

201 . 6 

1 .01 *10~5 

-3.97»105 

2016. 

4 .9  *1 o" 

10 

-40.32 

2 

Exact 

1260.1 

6.30x10_5 

-2.48*106 

1  2601  . 

0. 

-252.03 

Mesh-1 

1 391 .9 

6.96*10'5 

-2.74*106 

1  3904. 

3.8*10' 

10 

-310.24 

Mesh-2 

1 292.9 

6. 46*10~5 

-2.55*106 

12917. 

4 .6*10' 

1 1 

-265.25 

Mesh -3 

1292.9 

6.46*10  ^ 

-2.55*106 

12917. 

4.6*10" 

1  1 

-258.46 

3 

Exact 

3534.1 

-  4 

1  .77*10 

-6 .96  *1 06 

35341  . 

0. 

-706.82 

Mesh-1 

4727 .9 

2.36*10"4 

-9.31 *106 

47120. 

1 .5*10' 

10 

-1260.49 

Mesh-2 

3790.1 

1 .90*10~4 

-7.46*106 

37809. 

5 .8  *1  o' 

1 1 

-818.98 

Mesh-3 

3790.1 

1 .90*10_l4 

-7.46*106 

37809. 

5 . 8  *  1 0~ 

1  1 

-757.10 

fE-1*107;  v-0;  t-0.1;  p-0.000254;  L - 1 0 ;  b-1 
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the  du/dh  sensitivity  is  still  less  accurate  (15.9  percent  error 
for  the  third  mode)  than  the  frequency  (7.2  percent  error). 

Mesh  3  is  different  from  Mesh  2  only  in  the  specification  of 
the  nodal  coordinate  sensitivities,  which  new  are  specified  so 
that  all  nodes  move  proportionally  when  the  length  changes.  Mesh 
3  produces  results  which  are  essentially  exact  for  Mode  1,  and 
reduces  the  error  in  du/dL  by  half  for  the  higher  modes.  For 
Mesh  3,  the  error  in  all  of  the  sensitivity  results  is  generally 
no  larger  than  the  frequency  error  for  each  case  considered. 

7.3.6  Twisted  Plate  Frequency  Sensitivity 

Figure  24  shows  a  twisted  cantilever  plate  which  has  been 
used  extensively  for  the  comparison  of  natural  frequency 
predictions.  We  wish  to  compare  natural  frequency  sensitiv¬ 
ities  obtained  with  the  present  analysis  to  those  derived  from 
finite  differencing.  For  the  case  considered  we  take  E=107, 
v=0 . 30 ,  p=0. 00026;  the  dimensions  are  length  a=3 ,  width  b**l,  and 
thickness  h=0.050. 

For  the  present  analysis,  we  employ  a  6x6  mesh  of  bilinear 
elements,  which  is  adequate  for  the  first  few  modes.  As  evidence 
of  this.  Table  16  summarizes  the  first  several  modes  predicted 
for  a  twist  angle  of  0=30°.  Frequencies  obtained  from  NASTRAN69 
using  a  mesh  of  128  TRIA2  elements  are  tabulated  as  well,  and  the 
two  solutions  are  in  reasonable  agreement. 

Table  17  lists  the  computed  natural  frequency  sensitivities 
for  modes  1-4,  for  a  twist  angle  of  32*.  The  finite  difference 
estimates  shown  for  the  sensitivities  have  been  obtained  from 
separate  eigenvalue  solutions  performed  for  twist  angles  of  31.9° 
and  32.1*.  The  agreement  of  the  predictions  is  reasonably  good, 
and  is  quite  accurate  where  the  corresponding  sensitivity  is 
large  in  magnitude. 
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Table  16. 


Frequency  Comparison  for  30° 


Twisted  Plate 


t 


Mode 

Type 

NASTRAN 

Present 

1 

1  -B 

3.42 

3.20 

2 

2-B 

19.10 

19  .08 

3 

1  -T 

26.04 

26.52 

4 

3-B 

60.15 

61 .62 

5 

1  -EB 

73.00 

74.19 

6 

2-T 

78.50 

82.93 

Normalized  frequencies  are  X  -  wa  /ph/D 
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Table  17.  Frequency  Sensitivities  for  32°  Twisted  Plate 


Mode 

Type 

cu(  31  .9°) 

u>(  32°) 

io(32.1°) 

3U/301 ,2 

Au>/A03 

1 

1-B 

921.16 

920.23 

919.29 

-467.3 

-535.1 

2 

2-B 

5354.91 

5346.21 

5337.53 

-4668.9 

-4979.0 

3 

1-T 

9006.48 

901 7.26 

9028.02 

4927.2 

6170.8 

4 

3-B 

17141  .7 

17119.2 

17096.8 

-1  2791 .5 

-1 2862.9 

Note  that  the  variable  0  is  defined  in  radians. 

2 

3w/30  is  computed  directly  by  the  sensitivity  solution. 
36w/A0  is  computed  by  differencing  values  at  31.9°  and  32.1°. 
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The  relatively  low  accuracy  of  the  sensitivities  in  torsion 
are  thought  to  be  an  artifact  of  the  bilinear  element  used  in  the 
frequency  calculations.  Since  the  element  is  integrated  with  a 
single  point,  the  twisted  undeformed  geometry  is  not  reflected  in 
the  stiffness  computation  (other  than  in  "nodal  offsets"  which 
occur  at  each  node,  and  which  are  taken  into  account) .  High 
accuracy  for  twisting  modes  (and  presumably  for  the  corresponding 
sensitivities)  therefore  requires  a  relatively  fine  mesh. 

7.4  PROBABILISTIC  ANALYSIS  EXAMPLES 

Probabilistic  solutions  obtained  with  the  methods  described 
herein  are  described  in  this  Section.  It  should  be  recognized 
that  the  finite  element  calculations  performed  in  the  probabil¬ 
istic  analyses  are  limited  to  the  basic  (deterministic)  solution 
and  the  sensitivity  analyses,  so  that  the  numerical  behavior 
reported  for  the  sensitivity  analyses  of  the  previous  Section  is 
typical  of  the  probabilistic  solutions  as  well.  The  additional 
steps  of  performing  variance  and  percentile  calculations  complete 
the  process. 

7.4.1  Forced  Vibration  of  a  Cantilever  Beam 

The  cantilever  beam  shown  in  Figure  22  (see  Sections  7.3.2 
and  7.3.5)  is  subjected  to  a  uniform  pressure  load  which  varies 
sinusoidally  in  time,  q=-0 . 01* sin (ut) .  The  finite  element  model 
used  is  the  same  as  Mesh  1  of  Section  7.3.5,  so  that  the  first 
resonant  frequency  is  at  <j=202.3  Hz.  We  consider  forcing  fre¬ 
quencies  in  the  range  200<<J<205 ,  to  determine  the  steady-state 
response  behavior  of  the  beam  near  its  first  mode. 

Figure  25  shows  the  amplitude  of  the  end  deflection  versus 
forcing  frequency.  Note  that  the  tip  displacement  and  the  forces 
are  in  phase  for  frequencies  lower  than  the  natural  frequency, 
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and  out  of  phase  after  crossing  the  resonance.  In  the  neighbor¬ 
hood  of  the  natural  frequency,  a  step  of  0.125  Hz.  has  been  used 
for  the  forcing  frequency  in  generating  the  results  of  Figure  25. 

Statistical  parameters  selected  for  this  analysis  are  the 

7  5  .  -4 

elastic  modulus  (E=lxio  ;  a_=lxio  ),  mass  density  (p=2.54xio 

-5  .  E 

a p= 1x10  ),  thickness  (t=0.10;  <7t=0.005),  and  beam  width  (b=l? 

ab=0.001).  Figure  26  contains  plots  of  amplitude  sensitivity  for 
each  of  these  parameters,  which  is  obtained  as  a  by-product  of 
the  probabilistic  solution. 

Figure  27  depicts  the  variance  of  the  tip  displacement 
amplitude  of  the  beam  versus  forcing  frequency.  The  upper  curve 
represents  the  total  variance,  and  reflects  the  probabilistic 
variation  of  all  four  parameters  (E,  p,  t,  b) .  The  remaining 
four  curves  show  the  contributions  to  this  total  from  individual 
'parameters.  The  relative  magnitudes  of  these  curves  depend  both 
upon  the  parameter  sensitivities  (Figure  26)  and  the  variances  in 
the  statistical  parameters.  In  this  case  the  thickness  variation 
is  the  most  pronounced  effect,  followed  by  the  density  variation. 

In  Figures  28  and  29,  we  show  the  probabilistic  solution  in 
terms  of  percentile  (confidence)  levels.  Figure  28  presents  this 
data  as  a  family  of  curves  for  discrete  confidence  levels.  The 
same  data  are  used  to  construct  a  continuous  surface  in  Figure 
29,  with  the  second  independent  variable  corresponding  to  the 
confidence  level.  Note  that  at  a  particular  frequency,  the  plot 
should  always  indicate  an  amplitude  which  increases  monotonically 
with  the  confidence  level. 

Figures  30  through  32  contain  moment  amplitude  results  for 
the  beam,  presented  in  forms  similar  to  the  displacement  solu¬ 
tions  above.  Bending  moment  values  are  obtained  from  the  element 
center  nearest  the  root  section  (x=l). 
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Sensitivity  of  Amplitude  to  Modulus 


Chntiimr  Beam 
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Figure  26.  Amplitude  Sensitivities  for  Cantilever  Beam. 


Sensitivity  of  Amplitude  to  Width 
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Figure  26.  Amplitude  Sensitivities  for  Cantilevered  Beam  (Concluded) . 
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TIP  DISPLACEMENT  AMPLITUDE 


Figure  27.  Displacement  Amplitude  Variance  versus  Frequency. 


Fiqure  28.  Tip  Displacement  versus  Frequency  and  Confidence  Level. 


TIP  DISPLACEMENT  AMPLITUDE 


Figure  29.  Displacement-Frequency-Confidence  Level  Surface. 


ROOT  SECTION  BENDING  MOMENT 


Figure  30.  Moment  Amplitude  Variance  versus  Frequency. 


ROOT  SECTION  BENDING  MOMENT 
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Figure  31.  Root  Moment  versus  Frequency  and  Confidence  Level. 


ROOT  SECTION  BENDING 
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Figure  32.  Moment-Frequency-Confidence  Level  Surface. 


7.4.2  Natural  Frequencies  of  a  Twisted  Blade 

This  problem  considers  a  2x2  inch  blade  with  45°  twist,  for 
which  experimental  results  have  been  collected  at  the  Air  Force 
Wright  Aero  Propusion  and  Power  Laboratory,  Wright  Patterson  Air 
Force  Base.  The  finite  element  model  of  a  single  blade  is  shown 
in  Figure  33.  In  the  experiments,  data  were  obtained  from  a 
twelve-bladed  disk  machined  from  flat  stock  and  twisted  to  the 
final  shape.  Blade-alone  frequencies  have  been  measured  by 
clamping  the  disk  near  the  blade  roots,  and  forcing  each  blade 
with  a  magnetic  exciter. 

The  blade  has  inner  and  outer  radii  of  4  and  6  inches,  and 
is  a  uniform  0.078  inch  in  thickness.  The  actual  part  contains 
1/4-inch  holes  at  either  edge  of  the  blade  roots;  in  the  model, 
we  have  simply  moved  the  root  nodes  inward  to  obtain  the  correct 
area  and  moment  of  inertia  there. 

The  first  two  statistical  parameters  used  in  this  example 

are  material  properties.  The  part  is  made  from  steel,  for  which 

we  take  E=29xio6  and  a  =87000.  A  density  of  0.000751  lb-sec2/in4 

is  assumed,  with  a  =3.7xio 

P 

The  remaining  parameters  have  to  do  with  the  twist  angle  of 
the  blade.  Let  £=x/S,  a  normalized  spanwise  coordinate.  For  the 
angle  of  twist  at  the  centerline  of  the  blade,  we  assume  that 

9/9max  =  +  c^2  (177> 

For  C=0  the  angle  of  twist  varies  linearly  with  the  spanwise 
coordinate.  The  distribution  of  the  twist  angle  along  the  blade 
span  is  shown  in  Figure  34  for  various  values  of  the  parameter  C. 
The  third  and  fourth  statistical  parameters  are  the  maximum  angle 
of  twist,  0  ,  and  the  twist  parameter  C.  The  derivatives  of  8 

with  respect  to  these  variables  follow  from  equation  (177),  so 
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Element  Model  of  45-Degree  Twist  Blade. 


Normalized  Distance 


that  the  necessary  coordinate  sensitivities  are  relatively  simple 
to  define  in  equation  form. 

The  nominal  value  of  0  is  ir/4  (45*).  The  standard  devia- 

max 

tion  used  for  0  „  has  been  computed  from  the  tolerance  of  ±0.005 

max 

inch  on  blade  tip  height  above  or  below  a  fixed  reference  plane, 
giving  an  angular  tolerance  of  0.0087267  (0.5°).  A  nominal  value 
of  C=-l . 00 ,  with  ac=0.001,  is  used  for  the  twist  parameter. 

The  first  three  frequencies  and  corresponding  standard 
deviations  computed  with  the  probabilistic  model  are  summarized 
in  Table  18.  The  experimental  results  for  the  lowest  (1-B)  mode 
are  suspect,  since  other  experiments  using  a  bladed  disk  of  the 
same  design  resulted  in  first  bending  frequencies  in  the  neigh¬ 
borhood  of  330  Hz.  It  is  apparent  that  some  details  of  the  root 
conditions  are  not  represented  perfectly  in  the  model.  The 
remaining  modes  are  more  sensitive  to  the  blade  twist  profile,  as 
shown  in  Figure  35. 

Figure  36  shows  the  variances  of  the  first  three  natural 
frequencies,  as  well  as  the  contribution  of  each  statistical 
parameter  to  the  total.  Note  that  the  frequency  values  are 
listed  in  radians  per  second.  The  second  bending  mode  is  quite 
sensitive  to  the  total  twist  angle,  while  (from  Figure  35)  the 
twist  profile  is  relatively  unimportant.  Conversely,  the  first 
torsion  mode  is  influenced  less  by  the  total  angle  of  twist  than 
by  the  twist  profile  as  detemined  by  parameter  C.  The  small 
variance  assigned  to  parameter  C  prevents  it  from  influencing  the 
total  variance  of  the  natural  frequency. 
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Table  18.  Natural  Frequencies  for  45°  Twisted  Plate 


Mode  Type  m±Au)  (Experimental)  m±3o  (Computed) 

1  1-B  622.6  ±  1.3  Hz.  552.0  ±  29.4  Hz. 

2  2-B  1932.  ±  6.  2266.7  ±  278.3 

3  1-T  3335.  ±  1 .  3333.2  ±  236.4 
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Figure  35.  Blade  Frequencies  as  Functions  of  Twist  Parameter 
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Figure  36.  Frequency  Variances  for  Twisted  Blade. 
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APPENDIX  A 

PROTEC  INPUT  DATA  DESCRIPTIONS 


The  computer  program  in  which  the  analysis  techniques  re¬ 
ported  herein  are  implemented  is  called  PROTEC  (Probabilistic 
Sesponse  Of  Jurbine  IJngine  Components)  .  PROTEC  is  written  in 
ANSI  FORTRAN  77,  and  has  been  executed  successfully  on  CDC  Cyber, 
CRAY  X/MP,  and  DEC  VAX  machines.  This  Appendix  summarizes  input 
requirements  for  PROTEC.  The  remaining  Appendices  of  this  report 
describe  PROTEC  file  output  (Appendices  B  and  C) ,  data  conversion 
between  PROTEC  and  the  PATRAN10  modeling  package  (Appendix  D) , 
and  plotting  of  probabilistic  data  using  DISSPLA11  (Appendix  E) . 

Input  to  PROTEC  is  arranged  in  a  series  of  input  "blocks" . 
Each  input  block  begins  with  a  header  line  identifying  the  block, 
followed  by  the  data,  and  ends  with  a  blank  line  signifying  the 
end  of  the  block.  Input  block  types  are: 


Block  Name 

Status 

BOUNDARY 

Optional 

COORDINATE 

Required 

DERIVATIVES 

Optional 

DIAGNOSTICS 

Optional 

ELEMENT 

Required 

FORCE 

Optional 

GRAVITY 

Optional 

LAMINATE 

Optional 

MATERIAL 

Required 

OPTION 

Optional 

PARAMETERS 

Optional 

PRESSURE 

Optional 

PROPERTY 

Required 

TITLE 

Required 

Data  Description _ 

Nodal  boundary  conditions 
Nodal  coordinates 
Coordinate  derivative  data 
for  sensitivity  analysis 
Diagnostic  output  selection 
Element  connections 
Nodal  forces,  moments,  and 
prescribed  displacements 
Self-weight  loading 
Laminate  section  definitions 
Material  properties 
Analysis  options 
Statistical  and  sensitivity 
parameter  definitions 
Element  surface  pressures 
Element  thicknesses,  areas 
Alphanumeric  problem  title 


Input  blocks  may  appear  in  any  order  on  the  input  file.  While 
data  within  a  block  is  highly  structured,  comments  and  extra 
lines  may  be  inserted  between  blocks. 


A  -  1 


Formats  for  individual  input  blocks  are  described  on  the 
following  several  pages.  The  description  of  each  block  includes 

□  Header:  Four-character  block  title  (see  above) ; 

o  Format:  Typical  record  format  and  a  sample  data  line; 

□  Variables:  Definitions  of  input  variable  names;  and 

a  Notes:  Rules,  hints,  or  clarifications. 


A  -  2 


I 


BOUNDARY 
Input  Block 

BOUNDARY  Input  Block:  Nodal  Displacement  Constraints 


Header :  BOON 

Format : 


Variables: 

IBEG  =  First  node  number  to  be  constrained. 

IEND  =  Last  node  in  a  series  of  nodes  to  be  constrained. 
INCR  =  Node  number  increment. 

ID1-ID6  =  List  of  nodal  degrees  of  freedom  to  be  fixed;  the 
numeric  values  1-6  refer  to  u,  v,  w,  8  ,  8  ,  8 , 
respectively.  y 

Notes: 

o  If  IBEG,  IEND,  and  INCR  are  all  present,  each  of  the 
nodes  IBEG,  IBEG+INCR,  IBEG+2*INCR,  ...,  IEND  are 
constrained. 

o  If  IEND  is  omitted,  the  default  is  IBEG  (single  node) . 
o  If  INCR  is  omitted,  a  default  of  INCR  =  1  is  assumed, 
o  Only  one  degree-of-f reedom  value  (IDn)  is  required. 


A  -  3 


COORDINATE 

Input  Block 


COORDINATE  Input  Block:  Nodal  Coordinate  Data 

Header:  COOR 

Format: 


5  10 

20  . 

30 

40 

1  NODEl  INCR 1 

XCOORD 1 

YCOORD 1 

Z COORD  1 

Example : 


\  12 1  |  10 . 28 1  -67 . 42875 1  1.257E-2| 


Variables: 


NODE 
INCR 
XCOORD 
YCOORD 
Z COORD 


Current  node  number. 

Increment  for  node  number  generation. 

Cartesian  coordinate  X  at  the  current  node. 

Cartesian  coordinate  Y  at  the  current  node. 

Cartesian  coordinate  Z  at  the  current  node. 


Notes: 

o  Valid  node  numbers  are  from  one  to  the  maximum  number  in 
the  model.  Intermediate  node  numbers  may  be  omitted,  but 
must  be  constrained. 

o  INCR  is  used  for  generating  a  series  of  nodes  along  a 

line  from  two  successive  lines  of  data.  For  example,  the 
input 


10 

2.50 

-5.0 

3.10 

20 

2 

3.00 

-5.0 

2.10 

would  generate  the  following  coordinate  data: 


Node 

X 

Y 

Z 

10 

2 . 50 

-5.0 

3.10 

12 

2.60 

-5.0 

2.90 

14 

2.70 

-5.0 

2.70 

16 

2.80 

-5.0 

2.50 

18 

2 . 90 

-5.0 

2.30 

20 

3.00 

-5  “ 

2.10 
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DERIVATIVES 

Input  Block 


DERIVATIVES  Input  Block:  Coordinate  Derivatives  for  Sensitivity 

Analysis 


Header:  DERI  <value> 

Format: 

[  node |  incr|  xperiv)  yderiv|  zderiv| 
Example: 


tmt 

Variables: 

<value> 


NODE 

INCR 

XDERIV 

YDERIV 

ZDERIV 


l-~ . s.t§l  q.q| . -  laHI 


Integer  value  in  header  line,  specifying  IDENT 
(parameter  i.d.,  see  PARA  input  block)  for  the 
sensitivity  parameter  being  defined. 

Current  node  number. 

Increment  for  node  number  generation. 
Derivative  of  Cartesian  coordinate  X  at  the 
current  node,  with  respect  to  this  parameter. 
Derivative  of  Cartesian  coordinate  Y  at  the 
current  node,  with  respect  to  this  parameter. 
Derivative  of  Cartesian  coordinate  Z  at  the 
current  node,  with  respect  to  this  parameter. 


Notes: 

o  This  block  is  used  to  define  a  geometric  parameter  for 
use  in  sensitivity  analysis  (that  is,  a  parameter  which 
controls  Lho  placement  of  nodes  in  the  model) .  A  DERI 
block  is  required  for  each  such  parameter;  multiple  DERI 
blocks  are  distinguished  from  one  another  by  the  <value> 
appearing  in  the  header  line. 

o  NODE  numbers  are  as  defined  in  the  COOR  input  block,  and 
INCR  is  used  to  generate  data  exactly  as  the  COOR  block. 

o  If  the  sensitivity  parameter  is  p,  then  XDERIV  =  dX/dp, 
the  derivative  of  coordinate  X. 

o  Nodes  for  which  all  derivatives  are  zero  for  the  current 
parameter  may  be  omitted. 
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DIAGNOSTICS 
Input  Block 


DIAGNOSTICS  Input  Block:  Selection  of  Diagnostic  Output  Options 


Header :  DIAG 

Format : 


,  5  10- _ 15. 

| IDSW1 | IDSW2  IDSW3 


IDSW4 


|  IDSWn| 


Example: 


j\ . j°r . "f i  r~" . i . i . i . . 'i  i 


Variables: 

IDSWi  =  Number  of  a  diagnostic  output  switch  to  be 
activated  during  the  present  analysis. 


Notes: 


o  Continue  input  on  additional  lines  until  all  selections 
have  been  made.  Each  input  line  may  contain  from  one  to 
sixteen  switch  values. 


o  Valid  diagnostic  output  options  are  as  follows: 


IDSWi  Description  of  Diagnostic  Output 


1 

Element 

2 

Element 

3 

Element 

4 

Element 

5 

Element 

6 

Element 

7 

Element 

8 

Element 

9 

Element 

10 

Element 

11 

Element 

12 

Element 

data  (nodes,  properties,  coordinates) 
stiffness  matrices 
mass  matrices 

harmonic  stiffness  matrices  (K-XM) 
stabilization  (artificial)  forces 
transformation  matrices 
local  coordinates 
shape  functions 
strain-displacement  matrices 
stress-strain  matrices 
local  displacements 
displacement  sensitivities 
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ELEMENT 

Input  Block 


ELEMENT  Input  Block:  Finite  Element  Connections  and  Properties 

Header:  ELEM 

Format: 


Variables: 

Mnemonic  for  element  type. 

Element  number  for  current  element. 

Material  number  for  current  element.  A  negative 
value  refers  to  a  laminate  number,  as  defined  in 
the  LAMI  input  block. 

Physical  property  set  number  for  this  element. 
Node  increment  for  element  generation. 

Element  increment  for  element  generation. 

Nodes  connected  to  the  current  element,  listed  in 
counterclockwise  order  around  the  boundary. 

Notes: 

o  At  present,  the  only  acceptable  element  type  mnemonic 
(ETYP)  is  "SHELL",  designating  the  bilinear,  24-D.O.F. 
Mindlin  plate/shell  element.- 

o  Valid  element  numbers  are  from  one  to  the  total  number  of 
elements  in  the  model. 

o  Elements  may  be  generated  in  any  pattern  which  involves 
equal  increments  in  all  node  numbers  N1-N4.  INGEN  and 
IEGEN  appear  on  the  second  input  line  of  a  pair,  and 
specify  node  and  element  number  increments,  respectively. 
For  instance,  the  data 


generates  the  element  data: 


Element 

N1 

N2 

to 

N4 

20 

10 

14 

16 

12 

22 

13 

17 

19 

15 

24 

16 

20 

22 

18 

26 

19 

23 

25 

21 

ETYP 

IDEL 

MATL 


I  PR 
INGEN 
IEGEN 
N1-N4 
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FORCE 

Input  Block 


FORCE  Input  Block:  Imposed  Nodal  Forces,  Moments,  Displacements, 

and  Rotations 


Header :  FORC 

Format : 

1 - 5 - If - 

I  node]  kode]  VALUEj  ICASEJ 

Example: 


1  2581  T1  275.  ■5.1  I  .  g 


Variables: 

NODE 

KODE 


VALUE 
I  CASE 


=  Node  at  which  load  or  displacement  is  specified. 
=  Code  for  type  and  direction  of  prescribed  value: 


II 

<n 

»• 

2  =  F  ; 

3  =  F  ; 

4  =  M  ; 

5  -  M  :  6  = 

X 

y 

z 

X 

y 

7  =  u  ; 

8  =  u  : 

9  =  u  ; 

<*> 

II 

o 

H 

11=  0  ;  12= 

X 

y 

z 

X 

y 

=  Value  of  prescribed  force,  moment,  displacement, 
or  rotation. 

=  Static  load  case  number. 


Notes: 

o  The  first  three  values  must  be  provided.  There  are  no 
default  values  for  NODE  or  KODE.  ICASE,  if  omitted,  is 
assumed  to  be  1. 

o  If  a  nodal  displacement  or  rotation  is  set  to  zero  in 
this  input  block,  the  effect  is  the  same  as  a  constraint 
specified  in  the  BOUNDARY  input  block. 
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GRAVITY 
Input  Block 


GRAVITY  Input  Block:  Self-Weight  Body  Force  on  Entire  Model 


Header:  GRAV 

Format: 


30 
G  Z  | 


Example: 


-0  --0 1 


~o7o  | 


-386. \ 


Variables: 

GX,GY,GZ  =  Cartesian  components  of  gravity  vector,  defining 
both  the  magnitude  and  direction  of  the  local 
gravitational  acceleration. 


Notes: 

o  The  gravitational  force  per  unit  volume  at  any  point  is 
determined  from  F  =  p (GXi  +  GYj  +  GZk) ,  in  which  p  is  the 
material  density  at  the  point. 

o  By  default,  gravity  loads  become  part  of  load  case  1. 


A 


9 


LAMINATE 

Input  Block 


LAMINATE  Input  Block:  Laminate  Definitions  for  Layered  Shells 


Header:  LAMI 

Format : 


(1)  Sizing  Data  (one  per  laminated : 

10 


_ 5  10 

|  LAM |  NLAy| 

(2)  Laver  Data  (one  per  layer) : 

|  MATL^  THICK^  ANGLE^ 


Example: 


3 

0.060 

0.0 

2 

0.500 

45.0 

_ 1 

0.060 

0.0 

Variables: 

LAM  =  Laminate  number. 

NLAY  =  Number  of  layers  in  current  laminate. 

MATL  =  Material  number  for  a  specific  layer. 

THICK  =  Layer  thickness. 

ANGLE  =  Angle  from  local  'x'  axis  of  an  element  to  the 
material  '1'  axis  (fiber  direction). 

Notes: 

o  Laminate  definitions  must  be  numbered  sequentially  and 
input  in  ascending  order. 

o  The  layers  of  a  laminate  are  numbered  from  bottom  (layer 
1)  to  top  (layer  NLAY). 

o  MATL  may  reference  either  an  isotropic  or  orthotropic 
material,  as  defined  in  the  MATERIAL  input  block. 

o  ANGLE  is  positive  counterclockwise  when  viewing  an 
element  from  the  top. 

o  ANGLE  is  measured  in  degrees. 
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MATERIAL 

Input  Block 


MATERIAL  Input  Block:  Material  Properties  Data 

Header:  MATE 

Format: 


(1)  For  isotropic  materials  (one  line/material) : 


.  5  10 

20 

30 

40 

50 

.1  -MAT  ,////_  j _ 

- eL_ 

XNU  1 

RHO 

SY 

(2)  For  orthotropic  materials  (two  lines/material) : 


5  10 _ Z0 _ 30 _ 40 _ 50 _ 60 _ 70 


mm i 

El 

E2 

G12 

G13 

G2  3  1 

i im 

mnm 

RHO 

Cl 

nmzn 

mnm 

r imm 

Examples: 


|  l|  |  1 . E7 |  0 . 3 |  2 . 5E-4 |  10000. ] 


2 

25.  E6 

1.E6 

L  0.25! 

0.5E6 

0. 5E6 

0.2E6 

9.2E-5 

IMrMiM 

Variables: 

MAT 

E 

XNU 

RHO 

SY 

El 

E2 

XNU  12 
G12 
G13 
G23 

Cl,  C2 
Notes: 


=  Material  number  for  current  material. 

=  Extensional  modulus. 

=  Poisson's  ratio. 

=  Mass  density. 

=  Yield  stress. 

=  Extensional  modulus  in  material  direction 
=  Extensional  modulus  in  material  direction 
=  Major  inplane  Poisson's  ratio. 

=  Shear  modulus  in  material  (1,2)  plane. 

=  Shear  modulus  in  material  (1,3)  plane. 

=  Shear  modulus  in  material  (2,3)  plane. 

=  Failure  stress  constants. 


o  Materials  may  be  entered  in  any  order,  but  should  be 
numbered  from  1  to  the  total  number  of  materials,  with 
few  gaps. 


o  The  relationship  E  =  2G(l+y)  is  assumed  for  isotropic 
materials . 


o  Mass  densities  must  be  entered  in  units  consistent  with 
force,  length,  and  time  units  used  elsewhere  in  input. 

o  Constants  Cl,  C2  are  currently  not  used. 


A 


11 


OPTION 

Input  Block 


OPTION  Input  Block:  Selection  of  Solution  Options 

Header:  OPTI 

Format: 

(Enter  keywords  and  values  as  described  below. 
All  input  in  this  block  may  be  in  free  format.) 

Examples : 


HARMONIC  ANALYSIS _ 

FREQUENCY  10.  20.2.  23.  24. 

STATIC _ 

SENSITIVITY  STATIC _ 


Valid  Options  and  Keywords: 

EIGENVALUE  ....  Selects  natural  frequency  solution 

FREQUENCY  <values>  Defines  forcing  frequencies  for  steady- 

state  harmonic  solution 

HARMONIC  .  Selects  steady-state  forced  harmonic 

vibration  solution 

LOADCASES  ....  Defines  number  of  static  loading  cases 

MODES  .  Requests  a  specified  number  of  natural 

frequencies  in  an  eigenvalue  analysis 

SENSITIVITY  <name>  Requests  sensitivity  analysis  following 

a  basic  solution,  to  determine  response 
derivatives 

SSITERATIONS  .  .  .  Defines  the  maximum  number  of  iteration 

cycles  for  eigen’  alue  solutions 

SSTOLERANCE  .  .  .  Defines  the  relative  accuracy  tolerance 

used  to  test  eigenvalue  convergence 

STATIC  .  Selects  linear  static  solution 

Notes: 

©  Linear  static  analysis  normally  requires  the  STATIC  and 
LOAD_CASES  options. 

o  Steadv-state  harmonic  analysis  normally  requires  the  use 
of  HARMONIC  and  FREQUENCY  options. 

g  Natural  frequency  analysis  normally  requires  the  use  of 
EIGENVALUE  and  MODES  options. 

o  Valid  names  for  the  SENSITIVITY  option  are:  STATIC, 
HARMONIC,  and  EIGENVALUE. 

o  The  first  four  characters  of  each  keyword  (shown  in  bold 
above)  must  be  present. 


L2 


OPTION 

Input  Block 
(Continued) 

Defaults : 

o  LOAD_CASFS=  1,  if  STATIC  option  is  specified. 

O  SSITERATIONS=  max(  2*MODES,  10  )  if  EIGENVALUE  specified. 
O  SSTOLERANCE=  l.E-6,  if  EIGENVALUE  specified. 
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PARAMETERS 
Input  Block 


PARAMETERS  Input  Block:  Definition  of  Control  Parameters  for 

Statistical  or  Sensitivity  Analysis 


Header :  PARA 

Format : 


(1)  Sizing  Data  (one  line  onlvl : 

|  npar|*~ 

(2)  Control  Parameter  Data  (one  line/pararoeter^ : 

_ 25 


_  10  15 

[  ipar| itype| IDENT| 


STDDEV 


Example: 


2 

1 

1 

4 

100000. 

2 

4 

999 

0.05 

Variables: 


NPAR  =  Number  of  control  parameters  to  be  defined. 

IPAR  =  Sequence  number  of  current  parameter. 

ITYPE  =  Parameter  type:  1  =  modulus;  2  =  density;  3  = 

thickness;  4  =  geometric. 

IDENT  =  I . D.  of  material,  property  set,  or  other  data 
corresponding  to  the  current  parameter. 

STDDEV  =  Standard  deviation  of  current  parameter. 

Notes : 

o  Valid  sequence  numbers  IPAR  are  from  1  to  NPAR;  numbers 
outside  this  range  are  ignored. 

o  For  ITYPE  =  1,2,3,  the  parameter  being  defined  is  simply 
a  property  value  defined  elsewhere  in  the  MATERIAL  data 
or  PROPERTY  data.  When  ITYPE  =  4,  the  parameter  controls 
the  positions  of  nodes  in  the  model,  and  requires  some 
additional  data  for  its  definition  (see  DERI  block) . 

o  IDENT  refers  to  a  material  number  if  ITYPE  =  1  or  2 .  If 
ITYPE  =  3,  IDENT  refers  to  a  physical  property  number  as 
defined  in  the  PROPERTY  input  block. 


PARAMETERS 

Input  Block 
(Continued) 


o  When  ITYPE  =  4,  the  geometric  parameter  is  defined  by  the 
derivatives  dx/dp,  dY/dp,  dZ/dp  of  coordinates  at  certain 
nodes.  These  derivatives  must  be  specified  in  a  DERIVAT¬ 
IVE  input  block,  with  the  value  of  IDENT  specified  in  the 
block  header. 

o  STDDEV  is  unnecessary  for  sensitivity  analysis  alone,  but 
must  be  defined  when  probabilistic  information  about  the 
response  is  to  be  computed. 

o  The  units  of  STDDEV  must  be  the  same  as  those  of  the  mean 
values  defined  elsewhere  (e.g.,  a  modulus  value  defined 
in  MATERIAL  data).  For  ITYPE  =  4,  STDDEV  might  have  the 
same  units  as  the  nodal  coordinate  data  (if  the  parameter 
is  a  key  dimension),  or  different  units  (if  the  parameter 
is  an  angle,  for  instance) . 
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PRESSURE 
Input  Block 


PRESSURE  Input  Block:  Element  Pressure  Loading 

Header:  PRES 

Format : 


10 


I  IEBEG [  IEEND| IENCR 


press]*  icase  [ 


Example: 


j _ 5.]_  . 3  5 1  _  HI  .  -?50.0j.  .  I] 


Variables: 

IEBEG  =  First  element  number  to  which  the  specified 
pressure  is  to  be  applied. 

IEEND  =  Last  element  to  which  pressure  is  applied. 
IENCR  =  Element  number  increment. 

PRESS  =  Surface  pressure,  positive  outward 

ICASE  =  Static  loading  condition  number. 


Notes: 

o  Pressures  are  applied  to  elements  IEBEG,  IEBEG+ IENCR, 
IEBEG+2* IENCR,  . . . ,  IEEND. 

o  If  IEEND  is  not  given,  its  default  is  IEBEG  (one  element 
loaded) . 

o  If  IENCR  is  not  specified,  the  increment  is  set  to  one 
(all  elements  from  IEBEG  to  IEEND  loaded) . 

o  The  "outward"  direction  for  an  element  is  determined  by 
the  ordering  of  its  nodes.  When  the  element  is  viewed 
from  the  top  (nodes  N1-N4  arranged  counterclockwise) ,  a 
positive  (outward)  pressure  acts  upward,  toward  the 
viewer. 

o  If  ICASE  is  omitted,  load  case  1  is  assumed. 
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PROPERTY 
Input  Block 


PROPERTY  Input  Block:  Element  Thicknesses  and  Areas 

Header:  PROP 

Format: 

_ 5 _ 15 

|  I PR |  VALUE) 

Example: 

|  4 |  0 . 375 | 

Variables: 

IPR  =  Property  set  number. 

VALUE  =  Property  value  (area  for  1-D  elements,  thickness 

for  2-D  elements  and  shells) . 

Notes : 

o  Property  sets  may  be  entered  in  any  order,  but  should  be 
numbered  from  1  to  the  total  number  of  distinct  element 
properties  (or  with  few  gaps) . 
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TITLE 

Input  Block 


TITLE  Input  Block: 

Header :  TITL 

Format: 

1 _ 80 

[title  .  .  _ __l 

Example: 

|  Sensitivity  Analysis  of  Blade  with  Variable  Twist 


Variables: 

TITLE  =  Alphanumeric  problem  title. 

Notes: 

o  TITLE  may  include  any  valid  alphanumeric  characters. 
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APPENDIX  B 

POSFIL  Results  Pile  Description 


This  Appendix  documents  the  results  file  output  written  from 
PROTEC.  The  results  file  POSFIL  is  a  formatted,  card-image  file 
whose  structure  is  rigid  (and  therefore  simple  to  read  from  other 
programs) .  The  PATRAN  translator  PROPAT  (see  Appendix  D)  is  an 
example  of  a  program  which  reads  this  results  file  and  transmits 
data  to  other  programs  for  analysis  and  display. 


Data  on  POSFIL  are  arranged  in  blocks,  similar  in  concept  to 
the  input  data  blocks  (Appendix  A) .  Each  data  block  begins  with 
a  header  line  identifying  the  block,  followed  by  the  data,  and 
ends  with  an  empty  line  signifying  the  end  of  the  block.  Types 
of  data  blocks  generated  as  output  include: 


Block  Name  Description _ 

BOUN  Nodal  boundary  conditions 

CORD  Nodal  coordinates 

DISP  Nodal  displacement 

DSEN  Nodal  displacement  sensitivities 

ELEM  Element  connections 

ESEN  Eigenvalue  (frequency)  sensitivities 

FREQ  Harmonic  forcing  frequencies  or  system  natural 

frequencies 

LOAD  Nodal  forces  and  prescribed  displacements 

MATL  Material  properties 

MSEN  Mode  shape  sensitivity  coefficients 

PATR  Patran  neutral  file  title 

PVAR  Sensitivity  parameter  variances 

REAC  Nodal  force  reactions 

SSEN  Element  stress  sensitivities 

STRS  Element  stress  resultants 

TITL  Alphanumeric  problem  title 


Formats  for  the  individual  data  blocks  and  block  headers  are 
summarized  on  the  following  pages. 


B 


1 


AB 


' BOUN ' 
NUMDOF 


DOFKOD 


DESCRIPTIONS 


DISP 


ELEM 


NODE 
DISP (6) 


' DSEN ' 

I  CASE 
NUMNOD 
PARAM 


NODE 
DISP ( 6 ) 


' ELEM ' 


\wm  Hi  >  iji 


ELTYPE 
IELT 
MATLNO 
I  PROP 
CON  ( 4 


I 


Block  identifier 
Number  of  nodes 


Node  number 

Cartesian  coordinates  X. 


Block  identifier 
Load  case/mode  number 
Number  of  nodes 


Node  number 

Nodal  displacements  and 
rotations 


Block  identifier 
Lcod  case/mode  number 
Number  of  nodes 
Sensitivity  parameter  number 


FORMAT 


A8 ,  18 


80A1 


A8 ,  18 


18,  3E16.8 


A8 ,  218 


I 


Block  identifier 
Number  of  e 


Element  type 

Element  number 

Material  number 

Property  number 

List  of  connected  node  points 


Mode  number 

Sensitivity  parameter  number 
Natural  frequency 
Freouencv  sensitivit 


CO 

< 

318 

18,  8X , 
3E16 . 8 , / , 
16X. 3E16.8 

A8  , 

18 

> 

03 

718 

A8 , 

218 

218 

2E16 

/ 

.8 

B 


2 


BLOCK 

VARIABLE 

DESCRIPTIONS 

FORMAT 

FREQ 

' FREQ ' 
NUMMOD 

I  ANAL 

_ 

Block  identifier 

Number  of  frequencies 

'2'  =  Natural  frequencies 

'3'  =  Harmonic  forcina  frea's. 

A8 ,  218 

MM  lill  HHI 

IT  HI  M  1  MB  JM 

5E16.8 

LOAD 

'LOAD' 

'FORC' 

'GRAV' 

' PRES ' 

Block  identifier 
'FORC'  for  nodal  force 
'GRAV'  for  gravity  loading 
'PRES'  for  Dressure  loadina 

4A8 

NODE 

CODE ( 6 

FORCE 

Node  number 

'1'  =loading/ imposed  displ. 

'0'  =no  prescribed  values 
Applied  force  or  displ.  value 
for  each  direction 

18,  2X , 
6A1 , 

3E16.8,/, 
16X, 3E16 . 8 

MATL 

■ 

Block  identifier 

Number  of  materials 

A8 ,  18 

i 

ELMAT ( 9 ) 

Material  number 

Material  property  list 
_ (Ei  .  E?  ,  V  i  ■>  .  Gi  •>  .  Gi  i  .Gi  3.0  .  Ci  iCi  ) 

18,  8X , 
4E16.8, 

/  5E16.8 

1 

'MSEN' 

NUMMOD 

NUMPAR 

A8 ,  218 

MODE 

IPARAM 

INDEX 

COEFF 

Mode  number 

Sensitivity  parameter  number 
Coefficient  number 

Mode  shape  sensitivity  coeff¬ 
icient 

318,  E16.8 

PATR 

'PATR' 

Block  identifier 

A8 

DATE 

TIME 

PATVER 

Date  neutral  file  generated 

Time  neutral  file  generated 
Patran  version  number 

A12 ,  A8 , 
A10 

PVAR 

'PVAR' 

NUMPAR 

Block  identifier 

Number  of  sensitivitv  oaram's. 

A8 ,  18 

IPARAM 

ITYPE 

IDENT 

VAR 

Sensitivity  parameter  number 
'1'  =  material  modulus 
'2'  -  material  density 
'3'  =  thickness 

M'  =  geometric  parameter 
Material,  property,  or  DERIV 
block  i.d.  for  this  parameter 
Parameter  standard  deviation 

318,  E16.8 

VAR 


s 


TITL 


NODE 
FORC ( 6 ) 


IONS 


Block  identifier 

Applied  loading  case  number 

Number  of  nodes 


Node  number 

Nodal  reaction  force 

and  nodal  reaction  moments 


FORMAT 


A8 ,  218 


I  DEL 
ELTYPE 
I  CASE 
IANAL 


IP ARAM 


EPSS ( 8 ) 
SIGS (8 ) 


' STRS ' 
NCASE 
NUMELT 


IDEL 
ELTYPE 
I  CASE 
IANAL 


Element  number 
Element  type 
Load  case/mode  number 
' 4 '  =  Static  sensitivity 
'  5 '  =  Frequency  sensitivity 
'6'  =  Harmonic  sensitivity 
Sensitivity  oaramet 


Ge 

Ge 


lilUiilstJ1 


A8 ,  318 


I8,A8, 318 


'TITL' 


TITLE 


Block  identifier 

Number  of  load  cases  or  modes 

Number  of  elements 


Element  number 

Element  type 

Load  case/mode  number 

'1'  =  Static  solution 

'2'  =  Natural  frequency 

'3'  =  Steady  state  harmonic 


Generalized  strains 
Generalized  stresses 
von  Mises  stresses 


Block  identifie 


Alphanumeric  problem  title 


5E16.8 


A8 ,  218 


18 , A8 ,218 


5E16.8 


A8 


80A1 


B 


The  example  below  shows  the  POSFIL  output  for  a  very  simple 
finite  element  model.  For  larger  models,  the  nodes,  elements, 
degrees  of  freedom,  and  other  data  are  repeated  as  required  in 
the  same  formats. 


TITL 

Natural  frequencies  of  square  plate,  inplane  motions  only,  one  element 
HATl  1 

1  0 . 1 0OOOOOOE +08  0. 10000000E+08  0.25000000E+00  0.40000000E+07 

0.40000000E+07  0.40000000E+07  0.25900000E-03  0. 10000000E+06 

CORO  4 

1  0.00000000E+00  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

2  0.10000000E+01  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

3  0.10000000E+01  0.10000000E+02  O.OOOOOOOOE+OO 

4  O.OOOOOOOOE+OO  0. IOOOOOOOE+02  O.OOOOOOOOE+OO 

ELEH  1 

SHEL  11112  3  4 

BOUN  24 

111111011111001111101111 
FREO  2  2 

0.15426565E+10  0. 16475772E+12 

0ISP  1  4 

1  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

2  -0.25236442E-01  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

3  -0.25236442E-01  0.10000000E+01  O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

4  O.OOOOOOOOE+OO  0.10000000E+01  O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

D1SP  2  4 

1  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

2  0.10000000E+01  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

3  0. 10000000E+01  0.25236442E-01  O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

4  O.OOOOOOOOE+OO  0.25236442E-01  O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

STRS  2  1 

1  SHEl  1  2 

-0.25236442E-01  0.10000000E+00  0 . 69388939E ■ 1 7  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  -0. 12610265E+03  0. 499684 74E+05 

O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO  0. 10006329E+07  0.10006329E+07  0. 10006329E+07 

STRS  2  1 

1  SHEl  2  2 

0. 10000000E+01  0.25236442E-02  -0. 13010426E-17  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  0.53366982E+06  0. 13467928E+06 

-0.11368684E-12  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO  O.OOOOOOOOE+OO 

O.OOOOOOOOE+OO  0.961 3900 7E+07  0.96139007E+07  0.96139007E+07 
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APPENDIX  C 

LAYSTR  Layer  Stress  File  Description 


The  usual  stress  output  from  PROTEC  consists  of  reference 
surface  strains  and  curvatures,  as  well  as  force  and  moment 
resultants  at  the  center  of  each  element.  For  layered  elements, 
the  program  generates  much  more  detailed  stress  data  defining 
point  stress  distributions  throughout  the  element  thickness. 
However,  this  data  is  quite  lengthy  for  large  models,  and  often 
must  be  plotted  for  correct  interpretation. 

When  layered  elements  are  present  in  a  finite  element  model, 
SAFE  generates  a  separate  output  file  containing  detailed  layer 
stresses,  which  may  be  printed  or  read  as  input  for  graphical 
postprocessing.  The  name  of  this  lamina  stress  file  is  LAYSTR 
(LAYer  STResses) .  On  VAX  computers,  the  file  is  saved  automatic¬ 
ally  on  the  current  directory;  on  CDC  and  CRAY  systems,  LAYSTR  is 
a  local  file  which  must  be  saved  at  the  end  of  an  analysis  job. 

The  LAYSTR  file  is  a  formatted,  80-column  card  image  file, 
with  a  simple,  highly  structured  format.  For  each  element  and 
loading  case  (or  mode) ,  the  file  contains  the  following  data: 


Line 

1 


2 


3 


Format 

418 

218,  E16.9 

5E16.9 


1. 

2. 

3. 

4. 
1. 
2  . 

3  . 
1. 
2. 
3. 

4  . 

5. 


Data  Description _ 

IDEL  -  element  i.d.  number 
ICASE  -  load  case  or  mode  number 
LAMNO  -  laminate  i.d.  number 
NLAYER  -  number  of  layers 
LAYER  -  current  layer  number 
MATL  -  material  i.d.  for  layer 
Z  -  thickness  coordinate 

SXX  -  stress  component  a 

SYY  -  stress  component  o 

SXY  -  stress  component 

SXZ  -  stress  component  axy 

SYZ  -  stress  component  axz 


The  following  points  should  be  noted  concerning  the  data  items 
described  above: 


C  -  1 


o  Elements  appear  in  sequential  order  on  the  file. 


o  For  each  element,  all  load  cases  or  modes  will  appear 
together,  in  ascending  order. 

o  Within  an  element  and  case,  layers  will  appear  sequen¬ 
tially,  in  decreasing  order  (top  to  bottom) . 

o  For  each  layer,  three  "Z"  stations  are  output,  since  the 
computed  transverse  shear  stresses  vary  parabolically 
within  each  layer. 

o  Stress  components  are  referred  to  the  element  local  axes 

A  segment  of  a  typical  LAYSTR  file  corresponding  to  a  single  ele 
ment  with  three  layers  is  listed  below. 


1113 
3  1  0.250000000 

-28358.0728  -28302.5974  269.531802  O.OOOOOOOOOE+OO  O.OOOOOOOOOE+OO 

3  1  0.237500000 


-26940.1691 

-26887.4675 

256.055212 

7.42952757 

8.23334692 

3 

-25522.2655 

1  0.225000000 
-25472.3377 

242.578622 

14.4780537 

16.0444709 

2 

0.178805642 

2  0.225000000 
-0.178156580 

0.315352208E-02 

14.4780537 

16.0444709 

2  2  O.OOOOOOOOOE+OO 

O.OOOOOOOOOE+OO  0.000000000E+00 

o 

o 

♦ 

UJ 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

14.4786154 

16.0450934 

2 

0.178805642 

2-0.225000000 

0.178156580 

0.315352208E-02 

14.4780537 

16.0444709 

1 

25522.2655 

1-0.225000000 

25472.3377 

-242.578622 

14.4780537 

16.0444709 

1 

26940.1691 

1-0.237500000 

26887.4675 

-256.055212 

7.42952757 

8.23334692 

1 

28358.0728 

1-0.250000000 

28302.5974 

-269.531802 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

m 

♦ 

o 

o 

O.OOOOOOOOOE 
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APPENDIX  D 

PATRAN  INTERFACES  (PATPRO/ PROPAT) 


This  Appendix  describes  the  data  translation  performed  by 
PATPRO  (PATRAN-to-PROTEC)  and  PROPAT  ( PROTEC-to- PATRAN ) .  PATPRO 
converts  a  finite  element  neutral  file  from  the  geometric  model¬ 
ing  program  PATRAN  into  a  standard  input  file  for  finite  element 
analysis  by  PROTEC.  PROPAT  transforms  a  PROTEC  results  file  into 
a  PATRAN  results  file  for  postprocessing.  PATRAN10  is  a  product 
of  PDA  Engineering  in  Santa  Ana,  California. 

The  modeling-analysis-postprocessing  cycle  begins  in  PATRAN, 
where  the  finite  element  model  is  generated.  The  completed  model 
is  written  (by  PATRAN)  to  a  PATRAN  Neutral  File.  A  Neutral  File 
is  a  card-image  text  file  which  contains  geometric  data,  node  and 
element  definitions,  properties  data,  loads,  constraints,  and 
model  identification  parameters.  From  the  Neutral  File,  PATPRO 
generates  most  of  the  data  required  to  perform  a  finite  element 
analysis  with  PROTEC. 

When  the  analysis  is  complete,  the  results  file  POSFIL  (see 
Appendix  B)  may  be  processed  using  PROPAT  to  create  plotting  data 
files  compatible  with  PATRAN.  Results  files,  together  with  the 
original  PATRAN  Neutral  File,  are  then  used  within  PATRAN  for  the 
graphical  display  of  stress  and  displacement  results. 

Both  PATPRO  and  PROPAT  are  written  in  ANSI  FORTRAN-77,  and 
are  operational  on  the  DEC  VAX  under  VMS  and  CDC  Cyber  under  NOS. 
Important  features  and  limitations  for  each  of  the  programs  are 
noted  in  the  paragraphs  below, 

PATPRO  ( PATRAN-to-PROTEC ) :  PATPRO  uses  the  PATRAN  Neutral 
File  to  generate  most  of  the  PROTEC  input  needed  for  an  analysis. 
PATRAN  data  types  which  can  be  translated  are  shown  in  the  Table 
below. 
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PATRAN 

Packet 

PROTEC 
Data  Block 

25 

TITL 

26 

PATR 

1 

COOR 

2 

ELEM 

3 

MATE 

4 

PROP 

6 

PRES 

7 

FORC 

8 

FORC 

Description _ 

Problem  title 
Model  identification 
Nodal  coordinates 
Element  connections 
Material  properties 
Physical  properties 
Pressure  loads 
Nodal  forces 
Nodal  displacements 


Notes  and 
Restrictions 


QUAD/4  (SHELL)  only 
Isotropic  materials 
Element  thicknesses 
Element  avg.  only 


All  nodes  present  in  the  PATRAN  model  are  translated  into 
PROTEC  format,  without  resequencing.  The  model  should  be  fully 
equivalenced  (i.e,  duplicate  nodes  eliminated)  in  PATRAN  before 
writing  the  Neutral  File.  We  also  recommend  the  node  renumbering 
facilities  in  PATRAN,  v/hich  are  extremely  effective;  the  RMS 
WAVEFRONT  criterion  is  most  appropriate  when  the  analysis  is  to 
be  performed  using  PROTEC. 

When  data  blocks  other  than  those  listed  above  are  needed, 
these  must  be  entered  manually  using  a  text  editor.  Examples  are 
the  OPTIons,  SENSitivity,  and  LAMInate  input  blocks. 

PROPAT  ( PROTEC- to- PATRAN^ :  PROPAT  processes  the  results  file 
(POSFIL)  generated  by  PROTEC,  and  produces  PATRAN-compat ible 
files  containing  nodal  or  element  results  "columns".  The  PATRAN 
results  files  are  binary  files  and  cannot  be  listed  or  printed; 
PROPAT  will,  at  the  user's  option,  generate  formatted  versions  of 
the  results  files  for  printing.  For  postprocessing,  both  the 
binary  results  files  from  PROPAT  and  the  original  PATRAN  neutral 
file  must  be  supplied  to  PATRAN.  Postprocessing  options  include 
plots  of  deformed  geometry,  stress  or  displacement  contours,  and 
color-coded  plots  of  key  element  or  nodal  results  from  PROTEC. 
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The  listings  which  follow  demonstrate  the  operation  of  the 
PATRAN  interface  programs,  and  show  the  types  of  data  which  are 
generated  at  each  stage  of  the  process.  The  table  below  gives  a 
summary  of  the  sample  listings. 

Listing  Title _  Description _ 

D.l  PATRAN  Session  File  Keyboard  input  to  PATRAN 

D.2  PATRAN  Neutral  File  Model  as  output  from  PATRAN 

D.3  PATPRO  Execution  Change  PATRAN  data  to  PROTEC  format 

D.4  PROTEC  Input  Data  Final  PROTEC  input  file 

D.5  POSFIL  Results  File  Results  file  output  by  PROTEC 

D.6  PROPAT  Execution  Change  results  file  to  PATRAN  format 

D. 7  Element  Results  File  Element  results  as  used  in  PATRAN 

D.8  Nodal  Results  File  Nodal  results  as  used  in  PATRAN 

D.9  PATRAN  Session  Interactive  postprocessing 


D 
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Listing  D.l.  PATRAN  Session 


Listing  D.2.  PATRAN  Neutral  File. 
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Listing  D.2.  PATRAN  Neutral  File  (continued). 
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Id 


Listing  D.2.  PATRAN  Neutral  File  (concluded) . 


patpao 


PA  T  .  P  *  0 
PATMAN  TO  PAOTtC  TMAHSWTOA 
PAf_PMO  TMAMSLATKS  A  PATMAN  MUTMAL 
riu  INTO  A  PNOTtC  IIW  FIU 


pumk  chttm  nc  patmam  mutmal  fiumpi*  « 

7  MUTUAL 


t  I 


PUMK  KNTKM  TM  PAOTKC  INPUT  FIUNANK  ..« 


TNi  TITU  Or  PATMAN  MUTUAL  PIU  IS  .  .  . 
SIMPLIFIED  SLAOK/DISC  SKCTOM 
THIS  PATMAN  MUTUAL  PIU  UM  CHEATED  AT  M.40.S 


ON  S7/10/1S 

TTK  nu 


FINN  PATMAN  UKXSION 


OF  NOOKS 

_  or  KUNKNTS 

NUMKM  Or  NATKRIAl  PMOPKMTIKS  • 
-  Or  PMVSICAL  PMOPKMTIKS  • 


l.S 


«S 

3* 

1 

0 


PMOCKM1MQ  PIKASK  WAIT  .  .  . 

1  .MS  CP  SKCOMS  KNKCUT10N  TINK. 


Listing  D.3.  PATPF.O  Execution. 
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Listing  D.5.  POSFIL  Results  File  (concluded). 
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Listing  D.6.  PROPAT  Execution 


Listing  D.7.  Element  Results  File 
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Listing  D.8.  Nodal  Results  File. 


PDA/PATRAN-G  RELEASE  1 .5 

CUSTOMER  -  WRIGHT  PATTERSON  AIR  FORCE  BASE 

FOR  INFORMATION  ON  NEW  FEATURES  IN  RELEASE  1.5 
OBTAIN  A  PRINTOUT  OF  FILE  INF015 


PLEASE  INPUT  THE  DEVICE  NAME  (OR  "REPORT") : 

>4C14 

INPUT  GO,  SES,  HELP,  OR  PDA/PATRAN-G  EXECUTIVE  DIRECTIVE 

>GO 

PATRAN  DATA  FILE?  KNEW  2. OLD  3. LAST 

>1 

PREPARING  THE  DATA  BASE  SUB-SYSTEM 

222  PARTITIONS  TO  BE  INITIALIZED: 

220 

200 

180 

160 

140 

120 

100 

80 

60 

40 

20 

MODE?  1. GEOMETRY  MODEL  2. ANALYSIS  MODEL  3 . DISPLAY  4. NEUTRAL  SYS.  5. END 

>SET,PH1 ,OFF 

"PHI"  IS  NOW  OFF  (WAS  ON  ). 

MODE?  1. GEOMETRY  MODEL  2. ANALYSIS  MODEL  3. DISPLAY  4. NEUTRAL  SYS.  5. END 

>SET, LABE, OFF 

P  HAS  El  LABELS  -  F  F  F  F  PHASE2  LABELS  -  FFFFFFFF 

PHASE3  LABELS  -  T  T  T  LOAD  LABELS  -  T  T  T  T 


Listing  D.9.  PATRAN  Session. 
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MODE?  1. GEOMETRY  MODEL  2. ANALYSIS  MODEL  3-DISPLAY  4. NEUTRAL  SYS.  5. END 
>4 

NEUTRAL  FILE?  1 .CREATE  OUTPUT  2. INPUT  MODEL  3. POST-PROCESSING  4. END 
>2 

INPUT  NEUTRAL  FILE  NAME 
> NEUTRAL 

DO  YOU  WISH  TO  OFFSET  ANY  NEUTRAL  INPUT  IDS?  (Y/N) 

>N 

LAYERED  BEAM  WITH  EDGE  STIFFENERS 

SHALL  WE  PROCEED  WITH  THE  READING  OF  THIS  FILE?  (Y/N) 

>Y 

READING  NODE  RECORDS: 

100 

READING  ELEMENT  RECORDS: 

READING  MATERIAL  PROPERTY  RECORDS: 

READING  PHYSICAL  PROPERTY  RECORDS: 

READING  PRESSURE  RECORDS: 

READING  DISPLACEMENT  RECORDS: 

READING  GRID  RECORDS: 

READING  LINE  RECORDS: 

READING  PATCH  RECORDS: 

READING  HYPERPATCH  RECORDS: 

READING  DATA  RECORDS: 

READING  GFEG  RECORDS: 

READING  CFEG  RECORDS: 


Listing  D.9.  Continued 
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NUMBER  OF  ITEMS  READ  FROM  NEUTRAL  FILE: 

NUM  NODE,  ELEM,  MATL,  PROP,  CORD,  PRES,  FORC,  DISP,  DEFO,  TEMPN,  TEMPE 
45  32  1  1  0000000 

NUM  GRID,  LINE, PATCH,  HPAT.DLINE,  DPAT.DHPAT,  LIST,  DATA 
602000000 
NUM  GFEG,  CFEG 
2  2 

MODE?  1. GEOMETRY  MODEL  2. ANALYSIS  MODEL  3. DISPLAY  4. NEUTRAL  SYS.  5. END 


NEUTRAL  FILE?  1 .CREATE  OUTPUT  2. INPUT  MODEL  3. POST-PROCESSING  4. END 


POSTPROCESS?  1. DEFORMAT IONS  2 . ELEMENT  QUANTITIES  3. END 

v  * 


INPUT  THE  KIND  OF  ATTRIBUTE  YOU  WISH  TO  SEE,  EG:  ID;  MIDjPID; 

TEJf>;  PRES;  DISP.N;  STRAIN, N;  STRESS, N;  VON.N;  CXUMN.N;  LIGHT;  NORMAL 
>COL, 1 4 

INPUT  THE  NAME  OF  THE  ELEMENT  RESULTS  FILE: 

>  ELERES 

DATA  WIDTH  -  15 

FILE  TITLE  -  EXAMPLE  PROBLEM 
PROTEC  ANALYSIS 
EIGENVALUE  RESULTS 

DATA  VALUES  RANGE  FROM  .123E*05  TO  .467E*06 
ASSIGNMENT?  1.AUTO  2. MANUAL  3. SEMI-AUTO  4. USE  CURRENT  LEVELS  5. END 


POSTPROCESS?  1 .DEFORMATIONS  2. ELEMENT  QUANTITIES  3. END 
>SET,SPECT, 15.1 ,2.3,4,5,6,7,8,9.10,11,12,13,14,15 

POSTPROCESS?  1 .DEFORMATIONS  2. ELEMENT  QUANTITIES  3. END 
>RUN,CONT,COL,14 

INPUT  THE  RESULTS  FILE  NAME: 

>  ELERES 

AVERAGING  COLUJtJ  5  OF  ELEMENT  RESULTS  FILE  AT  NODES. 
DATA  WIDTH  -  15 

FILE  TITLE  -  EXAMPLE  PROBLEM 
PROTEC  ANALYSIS 
EIGENVALUE  RESULTS 

DATA  VALUES  RANGE  FROM  .123E«-05  TO  .467E*06 

Listinq  D.9.  Continued 
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ASSIOiMENT?  1  .AUTO  2. MANUAL  3. SEMI-AUTO  4. USE  CURRENT  LEVELS  5. END 
>1 

ASSIOJED  CONTOUR  VALUE  CODES  FOLLOW: 

A  .41 10E+06  B  .  3815E+06  C  .3520E+06 

D  .  3224E+06  E  .2929E+06  F  .2634E«-06 

G  .2338E+06  H  .2043E+06  I  .1748E+06 

J  . 1 452E+06  K  . 1 157E+06  L  .8617E+05 

M  .5664E+05  N  .2711E+05 

A  SINGLE  COLUMN  NODAL  FILE  CALLED  "PATNOD  "  HAS  BEEN  PRODUCED. 

POSTPROCESS?  1 . DEFORMATIONS  2. ELEMENT  QUANTITIES  3- END 
>RUN,HIDE,CONT 

BEGINNING  PHASE-II  HIDDEN  LINE  PLOT  OF  ACCURACY  LEVEL  .20 
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POSTPROCESS?  1. DEFORMATIONS  2. ELEMENT  QUANTITIES  3 -END 
>STOP 

RESTART  DATA  BEING  WRITTEN  ON  87/09/29 
PDA/PATRAN  COW>LETED 


Listing  D.9.  Concluded 


APPENDIX  E 

DISSPLA  INTERFACE  (PRODIS) 


PRODIS  (PROTEC-to-DISSPLA)  is  an  output  processor  for  PROTEC 
which  performs  two  primary  functions: 

□  probabilistic  (variance)  computations 

□  presentation  graphics  using  the  DISSPLA11  library 

PRODIS  uses  the  results  file  POSFIL  (Appendix  B)  to  generate  x-y 
plots,  surface  plots,  and  histograms.  Data  used  in  PRODIS  plots 
also  can  be  written  to  separate  files  for  use  in  other  programs. 
PRODIS  is  written  in  ANSI  FORTRAN-77,  and  is  operational  on  the 
DEC  VAX  under  VMS  and  CDC  Cyber  systems  under  NOS. 

PRODIS  generally  allows  plotting  of  any  quantity  versus 
another,  although  some  combinations  are  best  suited  for  specific 
types  of  analysis.  Quantities  which  can  be  selected  for  plotting 
include: 

□  displacement  components  at  a  specified  node 

□  displacement  magnitude  at  a  specified  node 

a  maximum  displacement  for  a  collection  of  nodes 

a  principal  moment  for  a  specified  element 

□  von  Mises  stress  for  a  specified  element 

a  maximum  moment  or  stress  for  a  collection  of  elements 

a  harmonic  forcing  frequency 

In  some  cases,  it  is  desirable  to  plot  only  one  of  the  above 
quantities  for  a  series  of  load  cases  (static  analysis)  or  modes 
(natural  frequency  analysis) .  PRODIS  will  generate  histograms 
for  such  cases,  which  permits  an  easy  comparison  of  effects  from 
different  analysis  cases.  Results  from  steady-state  harmonic 
analyses,  with  forcing  frequency  as  an  independent  variable,  are 
typically  presented  as  x-y  plots  or  3-D  surfaces. 

Two  modes  of  presentation  are  included  in  PRODIS  for  display 
of  probabilistic  data.  In  static  or  natural  frequency  analysis, 


E  -  1 


variance  data  for  nodal  or  element  results  can  be  displayed  in 
histogram  form.  The  histogram  shows  variances  in  the  requested 
quantity  for  each  individual  statistical  parameter,  and  for  all 
parameters  combined.  Recall  that,  for  any  result  r  which  depends 
on  the  statistical  parameters  p^,  the  total  variance  is  (see 

Section  4.3): 

Var[r]  =  £  (f£-)2  Var[p.] 

i=l  *1 

In  effect,  the  histogram  displays  each  term  in  this  series  as 
well  as  the  total,  for  each  of  a  series  of  loading  conditions  or 
vibration  modes.  This  type  of  plot  is  useful  for  determining 
which  statistical  parameters  contribute  most  to  the  uncertainty 
in  the  computed  result,  and  for  comparing  this  data  for  different 
modes  or  loading  conditions. 

The  second  mode  of  presentation  for  probabilistic  results  is 
most  often  used  in  steady-state  harmonic  analysis,  where  forcing 
frequency  is  nearly  always  an  independent  variable.  This  being 
the  case,  one  can  assemble  frequency  response  (i.e.,  amplitude 
versus  frequency)  results  for  the  deterministic  response,  or  for 
a  given  percentile  level  (confidence  level) .  Amplitudes  versus 
both  forcing  frequency  and  confidence  level  may  be  presented  as  a 
family  of  curves,  or  as  a  three-dimensional  surface.  Some  plots 
of  this  type  can  be  found  in  Section  7.4. 

One  practical  concern  is  the  time  and  cost  associated  with 
processing  of  results.  The  results  which  are  generated  by  the 
basic  solution,  sensitivity  analyses,  and  probabilistic  computa¬ 
tions  often  represent  a  substantial  amount  of  numerical  data.  We 
recommend  using  the  "searching"  options  (those  which  search  for  a 
maximum  value  within  a  specified  set  of  nodes  or  elements)  with 
some  care,  since  a  great  deal  of  calculation  mav  be  required. 
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In  principal,  PRODIS  can  produce  output  on  any  graphical 
device  which  is  supported  by  DISSPLA.  However,  the  program  has 
been  tested  only  for  a  relatively  small  subset  of  these  devices. 
At  present,  PRODIS  is  equipped  to  generate  graphical  output  on: 

a  Tektronix  4000  series  graphics  terminals 

□  Calcomp  1051  drum  plotter 

□  Hewlett-Packard  7470-A  pen  plotter 

The  addition  of  other  DISSPLA-supported  devices  is  quite  simple, 
involving  only  a  call  to  the  appropriate  device  nomination  sub¬ 
routine  within  the  DISSPLA  library. 

PRODIS  is  fully  interactive,  and  issues  relatively  simple 
prompts  for  all  keyboard  input.  The  listing  which  follows  shows 
a  short  session  with  PRODIS,  and  the  resulting  histogram  plots. 
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Listing  E.l.  PRODIS  Sample  Execution. 


sting  E.l.  PROD IS  Sample  Execution  (continued). 
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Listing  E.l.  PRODIS  Sample  Execution  (continued). 


Listing  E.l.  PRODIS  Sample  Execution  (continued). 
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Listing  E.l.  PRODIS  Sample  Execution  (continued). 


E 


9 


I 


s 

i 


E  -  10 


Listing  E.l.  PRODIS  Sample  Execution  (concluded). 


